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PROGRAM  ABSTRACT 


TIVAR  is  a highly  modular  computer  program  for  generating  predictions  of 
ensemble  mean  and  standard  deviations  for  a closed-loop  marjual  control 
system.  The  system  dynamics  are  assumed  linear,  but  can  have  arbitrary  time 
variations.  The  feedback  control  is  generated  by  the  optimal  control  model 
of  human  response,  extended  to  treat  time-varying  systems  with  both  random 
and  deterministic  input  disturbances. 

The  TIVAR  program  is  written  in  the  FORTRAN-IV-EXTENDED  computer 
programming  language  and  is  designed  for  efficient  batch  operation  on  a 
Control  Data  CDC-6600  computer.  Data  input  to  the  program  is  provided  on 
standard  punched  cards  and  output  is  generated  via  the  lineprinter. 

In  this  manual,  we  give  the  problem  formulation,  the  methods  by  which 
system  parameters  can  be  changed,  the  input  deck  setup,  a description  of  the 
TIVAR  subroutines,  and  a sample  problem  with  its  solution.  The  manual  also 
includes  a description  of  and  instructions  for  using  the  auxiliary  gain 
computation  program  GNPROP. 
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2.  PROBLEM  DESCRIPTION 


2. 1 Svstem/DisDlav  Dynamics 

The  system  being  controlled  is  described  by  the  state-space  equation: 


x(t)  = Ax(t)  + Bu(t)  + Ew(t)  + Fz(t) 


where: 


x(t)  = system  state,  NX  vector,  of  which  NX1  states  are 
associated  with  input  shaping 
u(t)  = control  input,  NU  vector 

z(t)  = non-random,  deterministic  mean  or  bias  inputs,  NZ  vector 
w(t)  = independent  zero  mean  white  Gaussian  noise  inputs, 

NW  vector,  with  covariance: 

E[wi(t)  Wj,(a)]  = W5(t)  6(t-a)  i=1,...,NW  (la) 

The  system  parameters,  which  can  be  time  varying,  are: 

A = NX  by  NX  state  matrix 
B = NX  by  NU  control  matrix 
E = NX  by  NW  noise  matrix 
F = NX  by  NZ  mean  input  matrix 

The  displayed  system  variables,  which  includes  displayed  elements  plus 
their  rates  of  change  for  man-machine  studies,  are  given  by: 

y(t)  = Cx(t)  + Du(t)  (2) 

where : 

y(t)  = displayed  information,  NY  vector 
C = NY  by  NX  display  matrix  for  x(t) 

D = NY  by  NU  display  matrix  for  u(t) 

2.2  Feedback  Control  Solution 

The  system  (1)  is  assumed  to  be  controlled  to  minimize  a quadratic  cost 
f unctional : 


i 


2 


i 


1 
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J = E[x  (T  ) 5 x(T  ) + 1 / ^(x'q  X + y'Q  y + u'Q  u + u”Q  li)  dt]  (3) 

f ftT'i'  x y u r 

o 

where  and  T^.  are  the  initial  and  final  times,  respectively,  and  Q^,  Q^, 
and  are  diagonal  weighting  matrices.  In  the  program,  these  quantities  are 
read  and  stored  as  NX,  NY,  NU,  and  NU  vectors,  respectively.  The  control  law 
that  minimizes  J is  given  by: 


u(t)  = -L  x(t) 
c ~ 

for  continuous  time  implementation,  or: 


(4a) 


u = -L^  X 
n d “n 


(4b) 


for  discrete  time  implementation,  wherein  u is  treated  as  peicewise  constant 
over  intervals  of  length  = computer  time  step.  The  quantity  x is  a 
"besfestiraate  of  the  augmented  state  x = [x,u]'. 


The  NU  by  (NX+NU)  gain  matrix  or  may  be  input  to  the  program. 

The  logic  in  TIVAR  will  convert  from  one  to  the  other  as  necessary.  Note 

that  in  general,  L and/or  L will  vary  with  time.  An  interface  program 
c d 

GNPROP,  which  is  described  in  the  Appendix  of  this  manual,  provides  a means 

for  generating  time-varying  gains  L by  backwards  integration  of  the  Riccati 

d 

equation.  Alternatively,  TIVAR  contains  an  option  for  computing  the  feedback 

gains  L in  a suboptimal  manner.  The  assumptions  are: 
d 


1. 

2. 

3. 


T^-*-”,  i.e., system  time  constants  are  small  relative  to  the 
time-to-go. 

The  terminal  weightings,  Q = 0. 

is  computed  by  solving  the  steadv-state  discrete 
Riccati  equation,  using  the  program's  present  stored  values 


for  A,  B,  C,  Q 


and  Q 


X y u r 

The  continuous  time  feedback  control  law  (4a)  can  be  re-written  as: 


u(t)  = ■ 
or,  equivalently: 


[L 


cl 


c2 


LQ(t)J 
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T u + u = -L  .x(t)  + V (t)  (5) 

n opt  u 

where  v^Ct)  is  a "motor-noise";  is  a "neuro-motor"  time  constant.  The 

form  of  Eq.(5)  provides  a compatability  with  earlier  work  in  man-machine 

systems,  dealing  primarily  with  steady-state  problems  where  Q^,  is  adjusted  to 

give  ~ typically.  If  TIVAR  is  used  to  compute  L^,  the 

equivalent  continuous  gains  Lq  are  obtained  and  printed  along  with  the 

equivalent  L and  T . 

opt  n 

2.3  Human  Limitations 

The  human  generates  x(t)  on  the  basis  of  the  delayed  and  noisy  perceived 
information , 

y .(t)  = N ry.(t-T)]  ♦ V (t-T)  i=1,...,NY  (6) 

pi  i 1 yi 

T = human's  time  delay 


V (t) 

y 

and  N ( . ) 
i 

N^(x) 


observation  or  sensor  noise 


is  the  non-linear  observation  threshold. 
X - a.  X > a. 


0 

X + an 


|xl  < a^ 
X < -a^ 


(7) 


In  TIVAR,  Ni  is  replaced  by  its  Random  Input  Describing  function  which 
depends  on  the  mean  and  variance  of  the  input  signal  and  the  threshold  a^^. 
The  sensor  noise  i=1,...,NY  is  a zero-mean,  white  Gaussian  noise  with 

covariance: 


E[v  (t)  v'  (a)] 
yi  yi 


VO  (t)  , 


i=1 NY 


(8) 


that  contains  both  an  additive  and  a ratioed  component,  viz, 


V°i(t)  = Vy^(t)  + TT  E[y^2(t)]  (9) 

The  NY  additive  noise  variances  Vy^^  and  the  NY  noise/signal  ratios  Py^^ 

(in  dB)  are  inputs  to  TIVAR.  The  quantity  f^  is  the  attentional  allocation 
to  displayed  variable  y^.  The  f^  are  constrained  by; 


Report  No . 3^64 


Bolt  Beranek  and  Newman  Inc. 


f.  > 0 
1 

NY 

I = f = "Workload",  or  total  attention 

i TOT 

The  neuro-motor  interface  portion  of  the  model  is  given  by  Equation  5. 

The  motor  noises  v (t) , i=1,...,NU  are  zero-mean,  white  Gaussian  noises  with 
ui 

covariance 

^(t-0)  (10) 

that  also  contains  an  additive  and  a ratioed  component, 

^ ^ui  (11) 

The  NU  additive  noise  variances  and  the  NU  noise/signal  ratios  (in 
dB)  are  inputs  to  TIVAR. 


2.4  Nominal  or  Equilibrium  Path 

In  the  previous  equations,  the  state  x(t)  represented  the  deviations  of 

the  system  motion  from  some  nominal  path,  Xnom(t) • Similarly,  y(t)  is  the 

deviations  from  some  nominal  Y_^„(t).  In  most  applications  X and  Y are 

nom  noni  noni 

zero,  i.e.,  X=0  is  the  equilibrium  path.  TIVAR  can  be  used  to  generate  a 

soecific  nominal  path  U (t),  X (t).  Y (t)  that  meets  an  NTF  terminal 

nom  nom  nom 

condition: 


H X(T^)  + H U(T„)  + c = 0 
X f u f 

for  the  deterministic  system: 

X(t)  X A X(t)  + B U(t)  ; [X(Tq),  U(T^)]  = given 

while  minimizing  the  cost  functional: 


(12) 


(13) 


■’non  --  1^'  »r  <’"> 

o 

The  system  equation  (13)  is  the  same  as  Equation  (1)  less  the  external 
disturbances  w(t)  and  z(t).  The  cost  functional  weighting  Qp  is  the  same  as 
in  Equation  (3). 
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The  terminal  conditions  (12)  may  be  written  more  compac  ^ as 


H X(Tf)  = 0 = [H^  I 1 c]  X(T^) 

U(Tf) 

1 

The  NTF  by  (NX+NU+1)  matrix  H is  an  input  to  TIVAR. 
the  output 


In  computing  and 


= C X + D U 
nom  nom  nom 


it  is  assumed  that  the  matrices  A,  B,  and  QR  are  constant  over  the  interval 

[Tq,T^].  If  at  some  time  t,  these  quantities  change,  it  is  necessary  to 

resolve  the  optimization  with  T =t  and  X (t)  = present  value  of  X .If 

o nom  nom 

this  is  not  resolved,  the  terminal  conditions  (12)  will  not  be  satisfied.  It 
is  assumed  that  A,  B,  QR  will  remain  constant  over  [t,T^]  when  computing 
X 

nom 


TIVAR  discretizes  the  continuous  time  system  equations  for  the  time 
propagation  of  the  ensemble  statistics.  The  discretization  interval  A is  a 
program  input.  The  ensemble  statistics  include  the  mean  and  standard 
deviations  of  any  state,  output,  or  control.  The  mean  component  = E[x(t)] 
results  from: 


1.  the  deterministic  input  z(t) 

2.  the  nominal  path  X _ (t) 

nom 

The  standard  deviation  about  the  mean,  is  a result  of: 

1.  the  random  input  w(t) 

2.  the  human  sensor  noise  v (t) 

y 

3.  the  human  motor  noise  v^(t) 


Report  No.  3464 


Bolt  Beranek  and  Newman  Inc. 


TIME  VARIATIONS 


In  executing  TIVAR,  any  of  the  input  parameters  can  be  changed  by  the 
user  at  any  time  t.  There  are  two  methods  for  changing  parameters.  The 
first  is  via  card  inputs;  the  second  is  via  a user  written  routine  called 
INTNEW. 

3. 1 System  Codes 

In  TIVAR,  a maximum  of  32  system  elements  may  be  changed  at  any  given 
time  step.  Each  element,  or  parameter,  is  identified  by  a unique 
alphanumeric  code  and/or  an  index  number.  Table  1 defines  the  system  codes 
used.  When  a given  parameter  I is  changed  at  time  t,  TIVAR  sets  a flag, 
IFLAG(I)  equal  to  1 for  one  time  step.  The  implication  of  the  parameter 
change  is  then  addressed  via  the  internal  logic  in  TIVAR.  If  there  has  been 
no  parameter  change,  the  internal  flags  remain  at  their  nominal  zero  value. 

3.2  Changing  System  Parameters 

Parameters  may  be  changed  at  time  t via  external  or  card  inputs.  The 
alphanumeric  code(s)  is  used  to  identify  the  parameter(s)  being  changed  at  a 
selected  time.  The  input  required  for  each  code  is  given  in  Section  4. 

Parameters  can  be  changed  periodically  via  an  internal  subroutine 
INTNEW.  The  user  must  supply  his  own  code  — in  FORTRAN  IV  — to  use  this 
option.  For  any  specified  code  index  number,  I,  the  subroutine  INTNEW(KEY) 
is  called  once  every  NDT(I)  time  steps  with  KEY=I.  The  32  values  of  NOT  are 
inputs  to  TIVAR.  When  INTNEW  is  called,  IFLAG(KEY)  is  set  to  1.  The  user 
must  supply  the  manner  in  which  the  parameters  are  to  be  changed.  If  no 
FuPTn'iN  cjde  is  supplied  to  update  the  parameters,  no  changes  are  made. 

Thus,  TIVAR  assumes  that  the  "new"  parameters  are  identical  to  the  previously 
existing  ones. 

The  ability  to  change  parameters  easily,  via  user-written  code,  enables 
the  user  to  study  problems  in  which: 
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1.  system  variables  are  changing  continuously,  i.e.  evc.  y NDT=k 
time  steps. 

2.  only  certain  elements  in  the  system  matrices  are  changing  with  time 

3.  system  variables  are  functions  of  the  system  ensemble  statistics, 
e.g.,  nonlinear  systems  that  are  being  linearized. 

4.  logic  decisions  determine  how  parameters  change. 
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Table  1:  PARAMETER  CODES  IN  TIVAR 


CODE 

INDEX 

DESCRIPTION 

A 

1 

System  A matrix,  Eq.  (1) 

B 

2 

System  B matrix,  Eq.  (1) 

C 

3 

Output  C matrix,  Eq.  (2) 

D 

4 

Output  D matrix,  Eq.  (2) 

E 

5 

Noise  matrix,  E,  Eq.  (1) 

QX 

6 

State  weightings  vector,  Eq.  (3) 

QY 

7 

Output  weightings  vector,  Eq.  (3) 

QU 

8 

Control  weightings  vector,  Eq.  (3) 

QR 

9 

Control  rate  weightings  vector,  Eq.  (3) 

TD 

10 

Human's  time  delay  x,  Eq.  (6) 

PN 

1 1 

Variance  of  random  inputs  w,  Eq.  (la) 

MNA 

12 

Additive  motor  noise  variances,  V . Ea. 

u 

MNR 

13 

Motor  noise  ratios,  p^,  Eq.  (11) 

SNA 

14 

Additive  sensor  noise  variances,  V^,  Eq. 

SNR 

15 

Sensor  noise  ratios,  Py,  Eq.  (9) 

TH 

16 

Observational  thresholds,  a^,  Eq.  (7) 

ATTN 

17 

Attentional  allocations,  f^,  Eq.  (8) 

CGAIN 

18 

Continuous  time  control  gains,  Eq.  (4a) 

DGAIN 

19 

Discretized  control  gains,  Eq.  (4b) 

PNTVL 

20 

Time  interval  for  program  printouts 

NEWG 

21 

Trigger  to  compute  Lq  in  program 

XMINC 

22 

Increment  mean  value  of  states 

XSINC 

23 

Increment  standard  deviation  of  states 

F 

24 

Deterministic  input  matrix,  F,  Eq.  (1) 

XNOM 

25 

Nominal  path  initial  condition,  Eq.  (13) 

TCR 

26 

Terminal  condition  matrix,  H,  Eq.  (15) 

I NT 

27 

Transfer  control  to  subroutine  INTNEW 

DUMMY 

28-32 

Dummy  codes,  inactive  at  present 
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4.  INPUT  DECK  SETUP 

The  input  deck  structure  for  TIVAR  is  discussed  along  with  the 
user-written  routines  INTNEW  and  FDOT. 

4. 1 Control  Cards 

There  are  six  major  control  cards  required  by  TIVAR. 

Card  1 - Title  Information 
Column  1 : blank 

Columns  2-80:  alphanumeric  title  information 

Card  2 - Dimension  Information,  715  Format 

Field  1 : NX  = number  of  system  states 
Field  2:  NX1=  number  of  noise  shaping  states 
Field  3:  NU  = number  of  control  inputs  < 4 

Field  4:  NW  = number  of  random  noise  sources 

Field  5:  NY  = number  of  displayed  outputs 

Field  6:  NZ  = number  of  deterministic  inputs 

Field  7:  NTF=  number  of  terminal  conditions  for 

nominal  path  trajectory 

The  size  restrictions  are: 

NX  + NU  < NDIM 
NY  < NDIM 

where  NDIM  is  the  array  size  in  the  DIMENSION  statements  (15). 

Card  3 - Time  Information,  4E10.0  Format 

Field  1:  DEL  = discrete  time  step  interval  (sec) 

Field  2:  TO  = initial  time  (sec) 

Field  3:  TEND  = terminal  time  (sec) 

Field  4:  TEXT  = time  extension  to  propagate  covariance 
TEXT  sec  past  TEND. 

Card  4 - Printing/Plotting  Control,  3(2011)  Format 

Field  1:  Cols.  1-20  Print/Plot  codes  for  states  1-NX 


- 10  - 
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Field  2:  Cols.  21-40  Print/Plot  codes  for  outputs  1-NY 
Field  3:  Cols.  41-60  Print/Plot  codes  for  controls  1-NU. 

Each  column  of  an  associated  field  corresponds  to  one  state,  output  or 

control.  A single  integer  governs  the  printing  or  plotting  of  the  results 
for  that  variable: 

0,  or  blank  = no  printing  or  plotting  of  the  variable 

1 = print  mean  and  standard  deviation  vs.  time 
3 = plot  X and  Y ± a vs.  time 

X 

2 = print  and  plot  ensemble  statistics 

A maximum  of  5 states,  5 outputs  and  5 controls  may  be  printed  on  wide  paper. 
Narrow  paper  will  accomodate  only  3 variables  per  states,  outputs  or 
controls. 

Cards  B-6  - NSTEP(I)  for  internal  time  breaks  3215  Format 

The  32  fields  are  associated  with  the  32  parameter  codes  in  TIVAR  on  a 
one-to-one  basis.  The  I-th  field  is  associated  with  code  I.  NSTEP(I)  is  the 
frequency  (number  of  time  steps)  at  which  subroutine  INTNEW  is  called 
internally  with  KEY=I,  starting  at  time  TO.  Calling  INTNEW  with  KEY=I  sets 
IFLAG(I)=1  for  one  time  step.  The  actual  parameter  values  must  be  changed  by 
user-written  code.  If  no  code  is  supplied,  the  associated  parameters  retain 
their  previous  values. 

4.2  External  Parameter  Cards 

The  remaining  input  cards  are  used  to  change  system  parameters  via 
external  read-in  at  specified  times.  The  deck  set-up  follows  a standard 
form. 

Time  Card  - Cols.  1-4  Alphanumeric  TIME 

Cols.  11-20  Time  of  external  break,  E10.0  Format 
Code  Card  - Cols.  1-5  One  of  the  alphanumeric  codes  given 
Table  1. 

Parameter  Cards  - The  new  parameter  values  required  by  the 


associated  code. 


d 
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The  sequence  of  Code  Card  followed  by  new  parameter  values  .'epeated  for 
all  items  that  the  user  wishes  to  change  at  the  given  time.  To  change 
parameters  at  another  time,  input  a new  time  card  with  the  time  of  the 
desired  change,  followed  by  a code  card,  parameter  cards,  code  card,  etc. 

When  using  external  (card)  updates,  the  following  rules  must  be  observed: 

1.  Time  breaks  must  occur  in  increasing  order. 

2.  The  code  cards  can  occur  in  any  order. 

3.  The  parameter  cards  must  immediately  follow  the  associated 
code  card. 

4.  Parameter  cards  must  be  input,  as  the  program  expects 
to  read  in  new  values. 

5.  The  last  card  in  the  deck  is  an  end  card,  containing  the 
alphanumeric  END  in  cols.  1-3. 

Table  2 gives  the  required  input  for  each  of  the  active  codes,  as  well  as  the 
initial  parameter  values  that  are  set  internally  by  the  program,  prior  to 
t=T0. 

Data  is  entered  on  the  cards  in  8E10.0  Format,  i.e.,  in  floating  point 
fields  of  10  columns  with  a maximum  of  8 fields  per  card.  The  numbers  may  be 
either  in  fixed-point  (decimal)  format  or  in  scientific  (exponential)  format 
with  the  exponent  right-justified  in  the  field.  Matrices  are  entered  one  row 
at  a time.  If  a row  contains  less  than  8 entries,  the  remaining  fields  on 
the  card  are  left  blank.  If  a row  contains  more  than  8 entries,  continue  on 
a second  card  for  that  row.  A new  row  always  begins  on  a new  card.  Vectors 
are  enterd  in  similar  8E10.0  format:  the  first  entry  in  the  first  field,  the 
second  entry  in  the  second  field,  etc. 

4.3  User  Written  Routines 

The  two  us'er  written  routines  are  INTNEW(KEY)  and  FD0T(K,T).  The  purpose 
of  INTNEW  has  been  discussed  earlier  as  regards  parameter  variations.  The 
function  FD0T(K,T)  is  used  to  provide  the  time  history  of  the  deterministic 
inputs  2^(t).  Thus  at  time  t,  and  for  input  K,  K=1,...,NZ: 
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TABLE  2:  CARD  DATA  INPUTS  AND  INITIALIZATION 


CODE 

A 

INDEX 

INPUT  DATA 

] 

;nitial  value 

1 

A^j  ; 1=1, ...,NX,  j=1,... 

,NX 

A=0 

B 

2 

Bj^j  i 1 = 1 , . . . , NX,  j = 1 , . . . 
C^j  ! 1= 1 , . . . , NY , J=1,... 

,NU 

B=0 

C 

3 

,NX 

C=0 

D 

4 

D^ j 1 1=1,..., NY , J=1,... 

,NU 

D=0 

E 

5 

E^ j ! 1=1,..., NX , j=1,... 

,NW 

E=0 

QX 

6 

; 1=1,..., NX 

Q =0 

X 

QY 

7 

Qyl  J 1= 1 1 ... » NY 

Qy  = 0 

QU 

8 

; 1=1,..., NU 

Q =0 
u 

QR 

9 

Qpi  ; 1=1,..., NU 

Q^=0 

TD 

10 

T 

T =0 

PN 

1 1 

WO  ; 1=1 N W 

w=o 

MNA 

12 

; 1=1 NU 

V =0 
u 

MNR 

13 

; 1=  1 , . . . , NU 

p =-«dB 
u 

SNA 

14 

; 1=1,..., NY 

V =0 

y 

SNR 

15 

; i=1 NY 

0 =-oodB 

TH 

16 

a^  ; 1=1, ...,NY 

a=0 

ATTN 

17 

; 1=1, ...,NY 

f=1 

CGAIN 

18 

(L^)^j  ; 1=1 NU,  j = 1, 

’ 1=1, ...|NU,  j=1. 

. . . ,NX+NU 

L =0 
c 

DGAIN 

19 

. . . ,NX+NU 

PNTVL 

20 

h = printout  Interval 

h=0 

NEWG 

21 

no  Input  required 

— 

XMINC 

22 

6X^  ; 1=1, . . . ,NX+NU 

5x=0 

XSINC 

23 

(6a|)i^  ; 1=1 , . . . ,NX+NU 

6o2=0 

F 

24 

F^J  ; 1=1 NX,  J = i,... 

,NZ 

F=0 

XNOM 

25 

(Xnom)l  ; i=1,....NX^NU 

X =0 
nom 

TCR 

26 

; 1=1 NTF,  j=1,. 

KEY  In  12  Format 

. ,NX+NU+1 

H=0 

INT 

27 

— 
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FD0T(K,T)  = Zj^(T) 

If  NZ=0,  i.e.,  no  deterministic  inputs,  FOOT  should  return  0.0.  The  user 
must  supply  his  own  code  for  FOOT. 

4.i|  Miscellaneous  Comments 

The  following  comments  are  pertinent  to  the  use  of  TIVAR. 

1.  If  they  occur  at  the  same  time  instant,  external  (i.e.,  card  input) 
updates  take  precedence  over  internal  updates. 

2.  Reading  control  gains  or  takes  precedence  over  an  internal 
computation  request  (NEWG). 


1 
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5.  TIVAR  SUBROUTINES 

The  various  subroutines  that  constitute  the  TIVAR  package  are  discussed 
briefly,  along  with  the  named  COMMON  blocks. 

5. 1 Subroutine  Descriptions 

The  following  subroutines  are  used  for  the  covariance  propagation 
program.  Their  functions  are  described  below  briefly. 

1.  Subroutine  TIVAR.  Provides  the  major  logic  control,  dimensioning 
declarations  and  initialization  for  the  entire  package.  Virtually 
no  computation  is  done  in  TIVAR.  TIVAR  reads  the  6 control  cards. 

2.  Subroutine  UPDATE.  Reads  input  code  cards  and  the  subsequent 
data  input  cards.  Calls  INTNEW  periodically  as  per  NSTEP. 

3.  Subroutine  INTNEW.  User  written  routine  for  internal  updates. 

4.  Subroutine  GPFBGN.  Obtains  discrete  feedback  gains  L.  from  1 

d c 

input,  or  via  Internal  computations.  Converts  L to  1 and 

c d 

vice-versa  as  necessary.  Outputs  equivalent  T^,  and  L 

N opt 

5.  Subroutine  COVAR.  Performs  a onetime-step  propagation  of  the 
vector  mean  and  matrix  covariance  equations,  using  the  current 
values  of  system  parameters.  Updates  matrix  one  time  step. 

This  routine  is  the  heart  of  the  entire  package. 

6.  Subroutine  FINAL.  Generates  the  nominal  trajectory  and 

Y , when  H>fO. 
nom 

7.  Subroutine  INFORM.  Stores  the  information  for  printing/plotting 
on  a disk  file  for  subsequent  output.  Computes  max  and  min 
scaling  for  the  plot  routine. 

8.  Subroutine  PRINTR.  Prints  and  plots  output  information  as 
requested  on  control  card  4. 
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9.  Function  XGAIN.  Computes  Random  Input  describing  action  ^ 

for  threshold  nonlinearity.  ' 

i 

10.  Function  FDOT.  User  written  routine  to  generate  deterministic  ^ 

inputs  z^(t) . j 

11.  Library  Routines.  In  addition  to  the  above  routines,  the  program  | , 

requires  many  of  the  subroutines  from  the  linear  system  library.  j' 

ii 

5 . 2 Common  Block  Usage  1 1 

Named  COMMON  blocks  are  used  to  store  most  data  arrays  and  to  pass  | j 

information  among  the  various  subroutines.  The  dimension  declarations  are  1 I 

given  in  the  primary  subprogram  TIVAR. 

1.  /COMMUN/  NCODES,  PINTVL , NSP,  ICODES(32),  IFLAG(32),  NSTEP(32), 

LPRNT(60),  XMAX(60),  XMIN(60) 

NCODES  = Number  of  possible  system  codes  in  Table  1 (32) 

PINTVL  = Printout  interval,  read  in  code  20  (PSP) 

ICODES  = Alphanumeric  code  identifiers 

IFLAG  = 0 or  1 flags  to  indicate  parameter  changes  j 

NSTEP  = frequencies  for  internal  updates,  control  cards  5-6  ‘ 

LPRNT  = print/plot  control  read  on  card  4 

XMAX,  XMIN  = max  and  min  values  for  dynamic  variables  (computed 
internally) 

2.  /INOU/  KIN,  KOUT,  KPUNCH,  KDISK 

Logical  unit  numbers  for  input/output  devices 

3.  /PL0T1/ 

Required  input  for  lineprint  plot  subroutine 


4.  /MAIN1/  NDIM,  NDIM1 , C0M1  /MAIN2/  COM2 

Common  blocks  required  for  library  subroutines.  NDIM  = dimension  ( 

of  program  square  arrays. 


r 


m 
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5.  /C0MP1/,  /C0MP2/,  /C0MP5/ 

Common  blocks  used  for  internal  computations  and  internal  storage 
of  temporary  matrices. 

6.  /INPTS/  BD,  SAV,  AD 

Discretized  system  variables 

A 5 

AD  = discrete  system  augmented  matrix  = e ° 

BD  = discrete  system  B matrix  r e ° da  B^ 

SAV  = storage  for  the  last  NU  columns  of  A 

7.  /INPTX/  N,  NX,  NXl,  NU,  AC(1) 

Continuous  system  state  parameters. 

AC  = augmented  system  matrix  (codes  1 and  2) 

N = NX  + NU. 

8.  /INPTY/  NY,  C(1)  /INPTW/  NW,  E(1)  /INPTWD/  NZ,  F(1) 

Input  information  for  display  outputs,  random  inputs  and 
deterministic  inputs,  respectively. 

9.  /TIMES/  TIME,  DEL,  TO,  TEND,  TEXT,  TD,  NPRED,  EA(1) 

TIME  = current  value  of  time,  t 

DEL  = discrete  time  step,  5,  (card  3) 

TO,  TEND,  TEXT  = initial,  final  and  time  extension  (card  3) 

TD  = time  delay  x,  NPRED  = [x/6] 

EA  z (A  )NPRED  = exp(A  x) 
d c 

10.  /WEIGHT/  QX(30),  QY(30),  QR(30),  PSS(1) 

QX,  QY,  QR  = (augmented)  state,  output  and  control  rate  weightings 
respectively. 

PSS  = steady-state  Riccati  equation  solution. 

11.  /RATIOS/  PU(30),  VU(30),  PY(30),  VY(30),  TH(30),  GTH(30) 

PU,  VU  = motor  noise  ratios,  additive  components 

PY,  VY  = observations  noise  ratios,  additive  components 

TH,  GTH  = observation  thresholds,  RIDF  gain 
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12.  /GAINBK/  GAINS(l) 

GAINS  = state  feedback  gains 

d 

13.  /NOISE/  VU0(30),  VY0(30),  W0(30),  SIGMA( 1) 

VUO  = total  motor  noise  covariance  at  time  t 

VYO  = total  observation  noise  covariance  at  time  t 

WO  = random  input  noise  variance 

SIGMA  = filter  Riccati  equation  solution  at  time  t 

14.  /INCRE/  XINC(30),  XCINC(30),  ATTN(30) 

XINC  = mean  increment  (Sc 
XCINC  = SD  increment  6a 

X 

ATTN  = attentional  allocation,  f 

1 

15.  /NOMNL/  NTF,  XNOM(30),  YN0M(30),  TCR(1) 

NTF  = number  of  terminal  conditions 

XNOM,  YNOM  = present  value  of  X , Y 

nora  nom 

TCR  = terminal  condition  matrix,  H,  code  26 

16.  /COVX/  YEAR (30),  YC0V( 1 ) 

XBAR  = present  value  of  augmented  x(t) 

XCOV  = present  state  covariance  matrix 

17.  /COVY/  YBAR(30),  YC0V( 1 ) 

YEAR  = present  value  of  mean  output,  y(t) 

YCOV  = present  output  covariance  matrix 

18.  /COVXH/  XHN(30),  XHCOV(I) 

XHN  = present  value  of  mean  state  estimate  x(t) 

XHCOV  = present  covariance  of  state  estimate 

19.  /COVEF/  ENF(30),  ENP(30),  EC0V( 1 ) 

ENF  = present  value  of  mean  filtering  error,  5^j,(t) 

ENP  = present  value  of  mean  prediction  error,  b (t) 

P 

ECOV  = filter  error  covariance  at  time  t 
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20.  /COVXHE/  DUM(30),  XHEC0V( 1 ) 

DUM  = dummy  vector 

XHECOV  = cross-covariance  between  e^  and  x at  time  t 
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6.  SAMPLE  PROBLEM 

A sample  problem  illustrating  many  of  the  features  of  the  TIVAR  program 
is  given  in  this  section.  A description  of  the  problem  is  presented  first, 
followed  by  a listing  of  the  user  written  subroutines,  the  input  data  deck, 
and  a listing  of  the  output. 

6. 1 Sample  Problem  Description 

This  problem  analyzes  an  AAA  tracking  task.  The  controller  tracks  the 
azimuth  angle  of  a target  which  is  executing  a level  fly-by.  The  key  element 
illustrated  by  this  problem  is  the  use  of  the  FDOT  function  to  generate  the 
time  history  of  the  deterministic  input,  i.e.  the  azimuth  trajectory  of  the 
target. 

The  controller  is  explicitly  presented  with  a display  of  the  azimuth 
sighting  error,  and  is  assumed  to  derive  the  corresponding  error  rate.  His 
task  is  to  minimize  this  error  by  controlling  a set  of  rate-aided 
second-order  sight  dynamics,  his  control  being  a hand-crank. 

The  target  being  tracked  is  executing  a constant  speed  straight  and 
level  fly-by  of  44  sec  duration.  The  range  of  the  target  at  crossover,  R^, 

is  3000  ft,  and  the  speed,  V,  is  733  ft/sec  producing  a maximum  azimuth 
angular  velocity  of  about  14  deg/sec  at  crossover.  Initially,  22  sec  before 
crossover,  the  target  is  16,126  ft  from  the  crossover  point,  its  azimuth 
angular  position  is  -79.46  deg,  and  its  azimuth  angular  velocity  is  0.4683 
deg/ sec . 

The  system  states  are  defined  as  follows: 

= target  azimuth  angular  position  (dq^rees) 

x^  s target  azimuth  angular  velocity  (deg/sec) 

x^  = sight  azimuth  angular  position  (degrees) 

Xj^  = sight  azimuth  angular  velocity  (deg/sec) 

X = integral  of  the  control  input 
5 


j 

i 


i 
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It  is  assumed  that  the  controller  employs  a "constant  velocity"  model  of 
the  target  position.  Consequently,  the  state  space  equations  for  the  first 
two  states  (the  deterministic  states)  are: 

3Cl(t)  = X2(t) 
x^Ct)  = z(t), 

where  z(t),  the  deterministic  input  is  the  azimuth  angular  acceleration  of 
the  target.  Thus,  z(t)  is  given  by: 

z(t)  = -J(V/I1^)2 , ,,6) 

(1  + (D(t)/Rc)‘^)^ 

where 

V = the  speed  of  the  target  = 733  ft/sec 

R^  = the  range  of  the  target  at  crossover  = 3000  ft 

D(t)  = the  distacne  of  the  target  from  the  crossover  point 

= D + Vt  = -16,126  + 733t 
o 

The  FOOT  function  computes  z(t)  according  to  Eq.  16. 

The  transfer  function  relating  the  sight  position  to  the  control  input 
is: 


x^(s)  _ 64(s+1) 

u(s)  s(s2+i2s+64) 


( 1 + 1 / s)-^ 


64 


s +12S+64 


(17) 


Consequently,  the  state  space  equations  for  the  last  three  states  (the 
controllable  states)  are: 


= -64x^(t)  - 12x^(t)  + 64x^(t)  + 64u(t)  (18) 

X5  = u(t) 

The  displayed  outputs  are  the  azimuth  sighting  error  and  error  rate  and 
are  given  by: 
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y,(t)  = x,(t)  - xj(t) 

Y2(t)  = x^Ct)  - Xij(t) 

Regarding  the  human's  inherent  liraitatons,  the  observation  noise  to 
signal  ratio  (SNR)  and  motor  noise  ratio  (MNR)  are  set  to  the  nominal  values 
of  -20dB  and  -25dB,  respectively.  The  perceptual  time  delay  (TD)  is  set  to 
0.20  sec.  Observational  thresholds  are  set  at  0.05  deg  for  y^  (corresponding 
to  1$  of  the  field  of  view  of  the  gunsight) , and  0.025  deg  for  y^ 
(corresponding  to  a nominal  differential  threshold  for  motion). 

The  control  gains,  CGAIN,  computed  by  another  program,  were  chosen  by 
setting  the  cost  on  error  to  unity,  and  adjusting  the  cost  on  control  rate  to 
produce  a neuro-motor  lag  of  0.1  sec. 

Finally,  since  the  mean  initial  states  are  far  from  zero,  the  human's 
estimator  requires  some  time  to  settle  down.  The  presence  of  the  human's 
time  delay  further  compounds  this  problem.  To  provide  an  initializing 
transition  period  while  the  human's  estimator  settles  down,  the  sample 
problem  is  started  at  t=-6.0  sec,  at  which  time  the  target  is  20,524  ft  from 
the  crossover  point,  its  azimuth  angular  position  is  -81.684  deg  and  its 
azimuth  angular  velocity  is  0.2928  deg/ sec.  During  this  initialization 
period,  the  human's  time  delay  is  set  to  zero,  and  printouts  are  suppressed. 
At  t=0,  however,  the  time  delay  is  set  to  0.20  sec,  and  the  printout  interval 
is  set  to  1 sec. 

The  states  are  incremented,  XINC,  so  that  the  initial  angular  position 
and  velocity  of  the  target  are  correct,  and  the  human's  state  estimates  are 
incremented,  XHINC,  so  that  the  initial  error  and  error  rate  are  zero. 


oo  oo  ooooo 
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6.2  User  Written  Subroutines  for  the  Sample  Problem 

TIVP  - PROBLEM  DEPENDENT  SUBPROGRAMS  FOR  TIVARI 
INCLUDES: 

1 - FUNCTION  FDOT 

2 - SUBROUTINE  INTNEW 

3 - SUBROUTINE  ADJNOM 


FUNCTION  FDOT(K,T) 

DATA 

1 XO,  YO,  VO  /-16126.0,  3000.0,  733.0/, 

A=X/Y0 

B=1.0+A»A 

FD0T=R«D 

RETURN 

END 


1 


MUST  BE  SUPPLIED  BY  USER  - PROBLEM  DEPENDENT 
COMMON 

/COMMUN/  Ng^DE^^^^INTVL , NSP,  ICODES(32),  IFLAG(32) 

IFLAG(KEY)=1 

RETURN 

END 


C 

C 


1 


E 

M 

Q 


OM 


OUTPUTS  FOR  NOMINAL  PATH 
PROBLEM  DEPENDENT 
COMMON 

TIME,  , 

YBAR(30) 

NTF,  XNOM(30),  YNOM(30),  TCR(15,16) 


/TIMES/ 

/COVY/ 

/NOMNL/ 

CONTINUE 

RETURN 

END 
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6.4  Output  Listing  for  the 

P-I-D  CONTROLLER 
2-Dec-76  14:28 


ample  Problem 


DYNAMICS  READ  FROM:  PIDTI.DYN 
NO.  OF  TOTAL  SYSTEM  STATES  5 
NO.  OF  NOISE  SHAPING  STATES  2 
NO.  OF  CONTROL  SYSTEM  INPUTS  1 
NO.  OF  RANDOM  NOISE  SOURCES  1 
NQ-  QF  DISPLAYED  OUTPUTS  2 
NO.  OF  DETERMINISTIC  INPUTS  1 
NO.  OF  TERMINAL  CONDITIONS  0 


^^!6tooo 

TERMINAL  TIME  = 44.000 

EXTENSION  = 0.000 


0.050 


INTERNAL  TIME  BREAKS  INDEX  CODE  NDT 


EXTERNAL  UPDATE  AT  TIME  -6.000 

F MATRIX: 

0. 

1.000E+00 

0. 

0. 

0. 


CODE  F 


EXTERNAL  UPDATE  AT  TIME 

PU  VECTOR: 

-2.500E+01 

EXTERNAL  UPDATE  AT  TIME 


-6.000 


-6.000 


PY  VECTOR: 
-2.000E+01 


-2.000E+01 


EXTERNAL  UPDATE  AT  TIME  -6.000 


VECTOR: 

;.000E-02 


2.500E-02 


EXTERNAL  UPDATE  AT  TIME  -6.000 
TD  =0. 

EXTERNAL  UPDATE  AT  TIME  -6.000 


CODE  MNR 


CODE  SNR 


CODE  TH 


CODE  TD 


CODE  CGAIN 


CGAIN  MATRIX" 

-1.539E+0i  -3.303E+00  6.359E+00  6.404E-01  9.034E+00 

1 .OOOE+01 


EXTERNAL  UPDATE  AT  TIME  -6.000  CODE  XMINC 
XMINC  VECTOR: 

-8.168E+01  2.928E-01  -8.168E+01  2.928E-01 


EXTERNAL  UPDATE  AT  TIME  -6.000  CODE  PNTVL 
PNTVL=  6.000E+00 

EQUIVALENT  DISCRETE  GAINS  GENERATED: 
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DGAIN  MATRIX: 
-1 .186E+01 
8.732E+00 


-2.867E+00 


4.083E+00 


4.556E-01 


J 

7.782E+00  ; 


EXTERNAL  UPDATE  AT  TIME  0.000  CODE  TD 

TD  = 2.000E-01 

EXTERNAL  UPDATE  AT  TIME  0.000  CODE  PNTVL 

PNTVL=  1.000E+00 
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I 

( 


TIME 

0.00 

1.00 

2.00 

3.00 

4.00 

5.00 

6.00 

7.00 

8.00 

5.00 

16.00 

11.00 
12.00 

11 

15.00 

16.00 


STATE 
MEAN  STD 


7.95E+01 
7.9PE+01 
7.84E+01 
7.78E+01 
7.72E+01 
7.65E+01 
7.57E+01 
7.47E+01 
7.37E+01 
7.25E+01 
7 . 1 2E+0 1 
6.96E+01 
6.78E+OI 

5.98E+OI 

-5.59E+01 


19.00 

20.00 

22’.  88 

23.00 

24. ( 


.00 


.00 

.00 


29.00 

36.00 

31.00 

32.00 


39.00 

40.00 

tm 

43.00 

44.00 


-3.66E+01 
■2.65E+01 
-I.43E+QI 
.6.23E-01 
1 .3IE+OI 
2.56E+OI 

UM] 

5.05E+01 

5.56E+OI 

6.29Etoi 

6.55E+01 

6.77E+01 

,:8] 
.26E+01 
7.38E+01 
7.48E+01 
7.57E+01 
7.65E+01 
7.73E+01 

7.91  E+01 
7.96E+01 


1 

DEV 
8.56E-O2 
8.53E-O2 
8.5OE-O2 
8.48E-02 
8.46E-02 
8.44E-02 
8.42E-02 
8.40E-02 
8.38E-02 
‘ 5E-02 

JE-02 
IE-02 
8.29E-02 

^i:8i 

8.23E-02 

8.20E-02 

hmi 

8. 14E-02 
8.12E-02 

?i:8i 

8.05E-02 

8.03E-02 

?:8^|:8I 

7.97E-02 
7.95E-02 
S-02 
■02 
-02 
7.87E-02 
7.86E-02 
7.PE-02 
7.82E-02 
7.80E-02 
7.79E-02 
7.77E-02 
7.76E-02 
7.74E-02 

7.70E-02 

7.69E-02 


STATE  2 STATE  3 


MEAN 

STD  DEV 

MEAN 

STD  DEV 

4.68E-01 

1 .06E-04 

-7.95E+01 

7.96E-O2 

5.14E-01 

5.64E-01 

1.06E-04 
1 .06E-04 

-7.9PE+01 

-7.PE+01 

7.93E-O2 

7.92E-02 

7.89E-02 

6.22E-01 

1 .06E-04 

-7.78E+OI 

6.89E-01 

1 .07E-04 

-7.72E+01 

7.85E-02 

7.67E-01 

6.60E-01 

1 .07E-04 

-7.65E+01 

7.80E-02 

1.07E-04 

-7.57E+01 

7.75E-02 

7.69E-02 

9.7OE-OI 

1 .07E-04 

-7.48E+01 

1 . 10E+00 

1.08E-04 

-7.37E+OI 

7.62E-02 

1 .26E+00 
1.46E+00 

1 .08E-04 
1.10E-04 

-7.26E+01 

-7.12E+01 

7.54E-02 

7.44E-02 

1 .7OE+OO 

1 . 10E-04 

-6.90E+OI 

7.3IE-O2 

2.00E+00 

1 .16E-04 

-6.78E+OI 

7.13E-02 

i;i?i:88 

:818?:8] 

8:KI;8i 

3.55E+00 

1.51E-04 

-5.99E+01 

5.66E-02 

4.42E+00 

1 .69E-04 

-5.6OE+OI 

4.15E-02 

^:^p:88 

9.06E+00 

l:PI:8il 

2.46E-04 

i5:iip:81 

-3.68E+01 

§:8fi:8l 

1 .03E-01 

1 . 1 2E+0 1 

3.29E-04 

-2.68E+01 

1 .33E-01 

1 .32E+QI 

1 .40E+01 

-1.45E+Q1 

-7.26E-OI 

1.'l^|l8l 

1 .33E+01 

4.53E-04 

1 .33E+OI 

1 .29E-02 

1.14E+01 

4.69E-04 

2.58E+OI 

8.46E-02 

5.65E+00 

i:8?i:8il 

5.70E-04 

j:tp:8l 

5.07E+01 

. 1.31E-01 

4.47E+00 

5.86E-04 

5.57E+01 

1.05E-01 

3.59E+00 

6.04E-04 

5.97E+01 

7.3IE-O2 

2.92E+00 

6.07E-04 

6.2^E+01 

3.7IE-O2 

2.41E+00 

6.08E-04 

6.5bE+01 

3.3IE-O2 

2.02E+00 

6.08E-04 

6.78E+OI 

5.O8E-O2 

1 .7IE+QQ 
1 .47E+00 

:88i:8H 

6.96E+QI 

7.12E+01 

5.88E-Q2 

6.32E-O2 

1 .27E+00 

6.09E-04 

7.26E+01 

6.57E-02 

1 . 1 1E+00 

6.09E-04 

7.38E+OI 

6.72E-02 

9.76E-01 

8.65E-01 

6.09E-04 

6.09E-04 

7.48E+01 

7.57E+01 

7.65E+01 

6.82E-02 

6.88E-02 

7.72E-01 

6.09E-04 

6.93E-02 

6.93E-01 

6.09E-04 

7.73E+01 

6.96E-02 

6.25E-01 

5.b?E-01 

t.mi 

6.98E-02 

6.99E-02 

5. 16E-01 

6.09E-04 

7. 91 E+01 

7.00E-02 

4.72E-01 

6.09E-04 

7.96E+OI 

7.0  IE-02 
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-7.955E+01  M 
-7-955E+01  S 
-7-955E+01  S 


n-Mm  r 

S-6.00E-I-00 


7.965E-t-01 
S 7.965E+01 
S 7.965E+01 


M 4.40E+01  + 
S 4.40E+01  + 
S 4.40E+01  +H 


STATE 


1 


28 


COCO 


SCO  CO 
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5i:§8 

-6.00E+00 


4.676E-01  M 
4.676E-01  S 
4.676E-01  S 


M 1.400E+01 
S 1.400E+01 
S 1.400E+01 


+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

•f 

+ 

+ 

+ 

+ 

+ 

+ 

s 

+ 

+ 

+ 

s 

+ 

+ 

+ 

s 

+ 

+s 

+ 

+s 

+ 

+ 

+ 

+s 

+ 

+ 

+ 

+ s 

+ 

+ 

+ 

+ s 

+ 

+ 

+ 

+ s 

+ 

+ 

+ 

+ s 

+ 

+ s 

+ 

+ 

+ 

+ s 

+ 

+ 

+ s 

+ 

+ 

+ 

+ s 

+ 

+ 

+ 

h++H 

h++ 

++++H 

^■+++++++^ 

K+++++++H 

h++++++-M* 

+ s 

+ 

+ 

+ 

+ 

s 

+ 

+. 

+ 

+ 

+ 

s 

+ 

+ 

+ 

+ 

s 

„ + 

+ 

S '*■ 

+ 

+ 

+ 

+ 

s . 

+ 

+ 

+ 

+ 

s + 

+ 

+ 

+ 

s 

s + 

+ 

+ 

+ 

s 

+ 

+ 

s+ 

+ 

+ 

s 

+ 

+ • 

s 

+ 

+ 

s 

+ 

+ 

+ 

+ s 

+ 

+ 

+ 

+ s 

+ 

+ s 

+ 

+ 

+ s 

+ 

+ 

+ 

+ s 

+ 

M 4.40E+01 
S 4.40E+01 
S 4.40E4-01 


-M-++J 

+ S 
+ s 
+ s 
+ s 
+s 
+s 
+s 
s 
s 


+ 

+ 

+ 

+ 

+ 

+ 


+ 

+ 

+ 

+ 

+ 

+ 


+ 


STATE 


29  - 


+ + + ++  + + ++ 


cows 
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:t:88i:88 

-6.00E+00 


M 4.40E+01 
S 4.40E+01 
S 4.40E+01 


STATE  3 


j 


i 


1 


I 


30 


coco 
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TIME 
0,00 
1 . 

2,( 

3.00 

4.00 


OUTPUT 
MEAN  STD 


2.47E-03 
5E- 
)E- 
J.20E-0: 
3.66E-03 


.00 

19.00 

20.00 
21.00 
22.00 
23.00 

24.00 

25.00 

26.00 

27.00 

26.00 

!3:88  - 

31.00 

32.00 


7.62E-0- 

9.9IE-03 

umi 

2.20E-02 

2.89E-02 
3.8IE-O2 
5.06E-02 
6.77E-02 
9.10E-02 
1 .22E-01 
1.62E-01 
2.06E-01 
2.34E-01 
2. 10E-01 
1 .03E-01 
1 .23E-01 
2.07E-01 

1 .77E-01 
1 .39E-01 


-5.89E-02 
-4.40E-02 
3E- 


•1.83E-02 
■1 .35E-02 


-4.38E-O: 
■2. 5 IE-03 

■VMtU 

1.12E-03 
1 .87E-03 


1 .53E-01 
1.32E-01 


OiE 

3.49E-02 

3.40E-02 


OUTPUT  2 
MEAN  STD  DEV 
-8.1  IE-03  6.29E-02 

1.72E-03  6.58E-02 

2.5IE-03  6.96E-O2 


4.99E-0- 

6.22E-03 


■95E-02 
1.23E-02 
4.62E-02 
5.I8E-O2 
5.97E-02 
7.O8E-O2 
8.60E-02 
1,06E-01 
1 .31E-01 
1 .56E-OI 
1.67E-01 
1 .40E-01 
7.95E-02 
1 . 17E-01 


1 .27E-0L 
1 .65E-02 
2.XIE~02 
2.8^E-02 
3.88E-02 
5.22E-02 
6.92E-02 
8.68E-02 
9.35E-02 

6.43E-O2 
-2.47E-02 
-1 .51E-01 
-2.26E-01 


1 .08E-01 

Um] 

4.36E-03 
1.1  IE-02 


mm 

7.16E-02 
6.0  IE-02 


4.24E-02 

3.96E-O2 


J.20E-02 

3.16E-02 


6.36E 


3.27E-04 
T.O^E-04 
-4.65E-05 
-1 .55E-04 

-3. 18E-04 
-3.7OE-O4 


3E-02 
9.90E-O2 

]:im] 

1 .45E-01 
1 .68E-01 
1 .97E-OI 
2.34E-OI 
2.82E-01 
3.47E-01 
4.32E-OI 
5.42E-01 
6.74E-01 
8.02E-01 
8.53E-OI 
6.60E-01 
4. 19E-01 
5.98E-OI 

m] 

7.73E-01 

6.54E-01 

mm 

3.29E-01 

2.65E-01 

ujm- 


1 

.53E-01 
•32E-01 
. 16E-Q1 
.03E-01 
.22E-02 
■37E-02 


1.57E-02 

6.15E-02 
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3.780E-01 

M 

N.  3 

3.780E-01 

S 

S 3 

3.78OE-OI 

S 

S 3 

J94E-01 

394E-01 


M-6.00E+00 
S-6.00E+00  + 
S-6.00E+00  + 

+ 


+ 

+ 


♦ 

+ 


+ 

+ 

+ 

+ 

+ 

+ 

S 

s 

•f 


+ 

+ 


-f 


+ 

•f 


M 


M 


M 4.40E+01  + 
S 4.40E+01  + 
S 4.40E+01 


+ 

+ 


+ 

+ 

+ 

+ 

•f 
+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

M 

+ S 
+ S 
+ 

+ 

M+ 

+M 

S M 

S+ 

+S 

+ S 
•+ 

+ 

+ 

+ 

+ 

+ 

•f 


s 

s 


M S 
M S 


M S 
M S 

s«  i 

SMS 

SMS 

SMS 

SMS 

SMS 

SMS 

M 


+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

•f 

+ 

+ 

+ 

+ 

+ 


+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


S 

s 


M 


M 


S 

s 

s 


s 

M + S 

M+ 

S + M 

S + 

S + 

3 + 

M + 

+ 

+ 

+ 

+ 

+ 

+ 

-► 

+ 

+ 


S 

M 


M 


S 


. M 
5 M 
S M 
S M 
S M 
S M 
S M 
SMS 
SMS 


S 

S 

S 

s 

s 

s 


•f 

+ 

+ 

+ 

+ 


+ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 


OUTPUT  1 


32  - 


+ ++  + + + + + + + + ++  00+++  + + 
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M-6.00E+00 
S-6.00E+00  + 


-8.774E-01  M 
-8.774E-01  S 
-8.774E-01  S 


S-6.00E+00 


+ 

+■ 

+ 

+ 

+ 

+ 

+ 

+ 

+ 

s 

+ s 

+ 

+ 

+ s 
+s 
+■ 

+ 

+ 

+ 

+ 

+^ 

+ 

+ 

+ 


M 8.667E-01 
S 8.667E-01 
S 8.667E-01 

+ 

4- 

4- 

+ 

4- 

4- 

+ 

+ 

4- 

+ 

4- 

4- 

+ 

4* 

4- 

+ 

SMS 

4- 

4- 

+ 

SMS 

4- 

+ 

+ 

SMS 

4- 

4- 

SMS 

4- 

4- 

+ 

SMS 

4- 

+ 

+ 

SMS 

4- 

4- 

+ 

SMS 

4- 

4- 

SMS 

4- 

4* 

+ 

SMS 

4- 

4* 

+ 

SMS 

4- 

+ 

SMS 

4- 

4- 

+ 

SMS 

4- 

4* 

+ 

SMS 

4- 

+ 

+ 

SMS 

+ 

4- 

+ s 

M 

+s 

4- 

S 

M 

+ 

s + 

S + 

M 

4- 

s + 

■f 

M 

4- 

s + 

+ 

M 

4- 

s + 

+ 

M 

4* 

s 

+ 

M 

4- 

s+ 

4* 

M 

+ 

s + 

+ M 

S 

4- 

4* 

4- 

M 

4- 

s + 

+ 

M 

4* 

s + 

4- 

M 

4- 

s + 

4* 

M 

+ 

s + 

4- 

M 

4- 

s + 

4- 

M 

4- 

s + 

4- 

M 

4- 

s + 

s+ 

M 

4-  S 

4* 

+s 

M 

s+ 

4- 

§ i:im]  t 

S 4.40E+01 


+ 

+ 

+ 

+ 


+ 


M 
M. 

M 
M 
M 
M 

M i 
M . 
S M , 

s M : 

SMS 


S 

S 


+ 

-► 

+ 

+ 

+ 

+ 

+ 


+ 

+ 

+ 


OUTPUT 
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TIME 

0.00 

1:88 

3.00 

4.00 

5.00 
6.00 

7.00 

5.00 

q.oo 

10.00 

11.00 

12.00 

1.00 

14.00 

15.00 

16.00 

17.00 

15.00 

19.00 
20.00 
21.00 
22.00 
23.00 

24.00 

25.00 

26.00 


CONTROL 
MEAN  STD 


.00 

.00 


51.00 

52.00 


35.00 

36.00 


39.00 

40.00 

41.00 

i| 

4 


2. 19E-01 

5.66E-01 
6.32E-01 
7.04E-01 
7.86E-01 
a.82E-01 
9.97E-01 
1 . 13E+00 
1 .30E+00 
1 .51E+00 
1 .76E+00 
.09E+00 
2.50E+00 
3.04E+00 
3.75E+00 
4.69E+00 
5.94E+00 
7.55E+00 
9.52E+00 
1 . 1 6E+0 1 
1 .31E+01 
1.35E+01 
1 .25E+01 


1 

DEV 
2.78E-02 


50E+00 
7.02E+00 
5.57E+00 


2.9OE+OO 

2.39E+00 

f:?8i:88 

1 .45E+00 
1 .26E+00 

8.56E-01 

7.64E-01 

6.86E-01 


'2E-C 
2.00E-02 
2.08E-02 
2.20E-02 
2.36E-Q2 
2.55E-02 
2.78E-02 
3.O6E-O2 
3. 4 IE-02 
3.84E-02 
4.38E-O2 
5.06E-02 
5.96E-O2 
7.15E-02 
S.75E-02 
1 .09E-01 
1 .37E-OI 
1 .69E-01 
1.99E-01 
2.05E-01 
1 .50E-01 
1 . 14E-01 
1.53E-01 
1.89E-01 
2.Q0E-01 
1.84E-01 
1 .55E-01 

hmi 

7.98E-02 

6.50E-02 

3.99E-O2 

3.52E-02 

umi 

2.6IE-02 

2.42E-02 

2.26E-02 


5.11E-01  1.91E-02 


stow 


oooooonooo  OOOOOOO  ooooo  ooooooo  o ooooo 
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7.  TIVAR  LISTING 
NO  TABS,  COLONS,  OF  FORMFEEDS 


TIVAR  - TIME  VARYING  MAN/MACHINE  ANALYZER 
INCLUDES 

1 - BLOCK  DATA  TIVDAT  - INITIALIZES  VARIOUS  COMMON  BLOCKS 

2 - TIVMAIN  - CALLS  SUBROUTINE  TIVAR 

3 - SUBROUTINE  TIVAR  - PRIMARY  SUBPROGRAM 


ALSO  REQUIRES  THE  FOLLOWING  SUBPROGRAM  FILES 


TIVCMP  - COMPUTATIONS  FOR  TIVAR 
INCLUDES 

1 - SUBROUTINE  GPFBGN 

2 - SUBROUTINE  TERMG 

3 - SUBROUTINE  COVAR 

4 - SUBROUTINE  FINAL 


GENERAL  PURPOSE  FBK  GAIN  COMPUTATION 
COMPUTES  GAINS  BASED  ON  TERMINAL  COST 
PROPAGATES  COVARIANCES 
COMPUTES  NOMINAL  PATH  FOR  TERMINAL 
CONDITION 


TIVP  - PROBLEM  DEPENDENT  SUBPROGRAMS  FOR  TIVARI 
INCLUDES 

1 - FUNCTION  FDOT  - SPECIFIES  DETERMINISTIC  INPUT 

2 - SUBROUTINE  INTNEW  - PERFORMS  INTERNAL  UPDATES 

3 - SUBROUTINE  ADJNOM  - ADJUSTS  OUTPUT  FOR  NOMINAL  PATH 


nVIO  - I/O  FOR  TIVAR 
INCLUDES 


SUBROUTINE 

2 - SUBROUTINE 

3 - SUBROUTINE 

4 - SUBROUTINE 

5 - SUBROUTINE 


INTLET 

UPDATE 

OUTLET 

INFORM 

PRINTR 


SPECIFY  INTERNAL  UPDATES 
PERFORM  EXTERNAL  UPDATES 
SPECIFY  VARIABLES  FOR  OUTPUT 
DO  OUTPUT  FOR  A SINGLE  TIME  STEP 
PRINT  ALL  OUTPUT  AT  THE  END  OF  A RUN 


KPLOT  - LINEPRINTER  PLOTTING  PACKAGE 
INCLUDES 


1 

2 


7 - 


SUBROUTINE  KPLOT 
SUBROUTINE  ADJUST 
SUBROUTINE  QINIT 
SUBROUTINE  KPLOTC 
SUBROUTINE  QPLOT 
SUBROUTINE  PUCE 
SUBROUTINE  QPRINT 


BLOCK  DATA  TIVDAT 

INITIALIZES  VARIOUS  COMMON  BLOCKS 

COMMON 

1 /COMMUN/  NCODES.  PINTVL,  NSP,  ICODES(32),  IFLAG(32), 

1 NSTEPC32),  LPRNT(60),  XMAX(60),  XMIN(60) 

2 /INOU/  KIN,  KOUT,  KPTR,  KPUNCH,  KDISK 

3 /PLOT1/  NV.  NH,  NCPW,  LW,  XL,  XH,  YL,  YH,  NXES,  NDIR,  1ST, 

3 NGLV,  NGLH,  BSYM.  GSYM,  PSYM,  ND1,  ND2,  NOUT 

R /FILES/  KKB,  NAMCOV,  NAMLPT,  NAMINT,  NAMDYN,  nAmTMP,  NAMBUF 

DATA 

1 NH,  NXES,  NDIR,  NGLV,  NGLH,  BSYM,  GSYM 

1 /101,  1,  10,  20,  20,  1H+,  1H+/, 

2 NOUT,  ND1,  NCODES,  NV 

2/6,  1,  32,  51/, 

3 KIN,  KOUT,  KPTR,  KPUNCH,  KDISK 

3 / 5!  6 6’,  7,  8/, 

4 (ICODES(I),  1=1,32) 

4 /1HA,  1HB,  1HC,  1HD,  1HE,  2HQX,  2HQY,  2HQU,  2HQR,  2HTD, 

4 2HPN,  3HMNA,  3HMNR,  3HSNA,  3HSNR,  2HTH,  4HATTN,  5HCGAIN, 
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4 

4 


5HDGAIN,  5HPNTVL,  4HNEWG,  5HXMINC.  5HXSINC,  1HF, 
SHTCR.IhINT,  3h6xT,  3HQYT,  3*5Hd6mMY/ 


END 


4HXN0M,  ^ 
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PROGRAM  TIVMAIN 

1 (INPUT,  OUTPUT,  TAPE5=INPUT,  TAPE6=0UTPUT, 

2 PUNCH,  DISK,  TAPE7=PUNCH,  TAPE8=DISK) 
MAIN  PROGRAM 

CALLS  SUBROUTINE  TIVAR 


CALL  TIVAR 
END 
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SUBROUTINE  TIVAR 
PRIMARY  SUBPROGRAM 


DIMENSION 

1 W0RK(15,15),  TITLE(8) 


COMMON 
1 /COMMUN/ 


/INOU/ 

/PL0T1/ 


4 /MAIN1/ 

5 /MAIN2/ 

6 /COMP1/ 

7 /COMP3/ 

6 /C0MP5/ 

9 /INPTS/ 

A /INPTX/ 

B /INPTY/ 

C /INPTW/ 

D /INPTWD/ 
E /TIMES/ 

F /WEIGHT/ 


G 

H 

I 

J 

K 

L 

M 

N 

0 

P 

Q 


NCODES.  PINTVL.  NSP,  ICODESn?),  IFLAG(32), 
NSTEP(52),  LPRftT(60),  XMAX(60),  XMIN(60) 

KIN,  KOUT,  KPTR.  KPUNCH.  KDISK 

NV,  NH,  NCPW,  LW,  XL,  XH,  YL,  YH,  NXES,  NDIR,  1ST, 
NGLV,  NGLH,  BSYM,  GSYM.  PSYM,  ND1,  ND2,  NOUT 
NDIM  NDIMl,  C0MU15,15) 

COM2(15,15) 

U(15, 15) 

V(15.15), 

EG(15, 15) 

BD(15,4),  SAV(15.4),  A(1S,15) 

N,  NX,  NX1,  NU,  AA(15,15) 

NY,  C(15,15) 

NW’  E(15!155 

NWf),  F(15,15)  , . , 

TIME,  MIL,  TO,  TEND,  TEXT,  TD,  NPRED,  EA(15,15) 
QX(30),  QY(30),  QR(30),  PSS(15,15) 


COMMON 

/TCOST/ 

/NOISE/ 

/RATIOS/ 

/GAINBK/ 

/INC RE/ 

/COVX/ 

/COVY/ 

/COVXHE/ 

/COVXH/ 

/COVEF/ 

/NOMNL/ 


QXT(30),  QYT(30), 


PINV(15,15) 

SIGMA(15,15) 

TH(30),  GTH(30) 


vuo(30),  VY0(30).  wo(30).  sigma 
PU(30),  vunol,  FY(30),  VY(30), 
GAINS(15,15)  , , , 

XMINC(30),  XSINC(30),  ATTN(30) 
XBAR(30),  XCOVn5,15) 

YBAR(30),  YC0V(15,15) 

DUM(30),  XHECOV(15,15) 

XHN(30),  XHCOV(15,15) 

ENF(30).  ENP(30).  ECOV(15,15) 

NTF,  XN0M(30),  YNOM(30),  TCR( 15,16) 


C SET  NDIM 

NDIM=15 
NDIM1=NDIM-»-1 

C ZERO  THE  VECTORS  AND  MATRICES 

1 NPRED=0 

TD=0.0 

REWIND  KDISK 

DO  10  1=1,30 

XMINCa)=6.0 

XSINC(I)=0.0 

QX(I)=0.0 

QY(I)=0.0 

QR(I)=0.0 

woa)=o.o 

TH(I)=Q.O 

ATTN(I)=1.0 

xbarU)=o.o 

YBAR(I)=0.0 

XHNajrO.O 

ENF(I)=0.0 

ENP(I)=0.0 

XNOM(I)  = 0.0 

YN0M(I)=0.0 

PUTl)=0.0 

VU(l)=0.0 
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10 


12 

C 

1120 

50 

1125 


C 

1160 

1180 

1 

2 

S 

5 

6 

7 

C 

1200 


1225 

1 

2 

2 

C 


59 

C 

1310 


PY(I)=0.0 

VY(I)=0.0 

CONTINUE 

DO  12  I=1,NDIM 

DO  12  Jsl.NDIM 

TCR(I,J)=0.0 

SIGMA(I,J)=0.0 

F(I,J)=0.0 

C(I,J)=0.0 

EU,J)  = 0.0 

AA(I,J)=0.0 

GAINS(I,J)=0.0 

XC0VTl,J)=0.0 

YC0V(1  J)=10000.0 

XHEC0V(I,J)=0.0 

XHCoy(i,j)=o.o 

EC0V(I. J)=1 .OE-16 

CONTINUE 

SPECIFY  THE  TITLE 

READ  (KIN.1120).  (TITLE(I),  1=1,8) 

IF  (e6f(kin))  46o,50 
FORMAT  (8A10) 

CALL  PAGEFD(KPTR, 1) 

WRITE  (KPTR, 1125;,  (TITLE(I),  1=1,8) 

FORMAT  (/,1H  ,8A10) 

CALL  DAYTIM(KPTR) 

CALL  LINEFD(KPTR,2) 

GET  THE  PROBLEM  DIMENSIONS 

READ  (KIN,1160)  NX,  NX1,  NU,  NW,  NY,  NWD,  NTF 
FORMAT  (1615) 

WRITE  (KPTR,1180)  NX,  NX1,  NU,  NW,  NY,  NWD,  NTF 
FORMAT  (28H  NO.  OF  TOTAL  SYSTEM  STATES  ,13,/, 
29H  NO.  OF  NOISE  SHAPING  STATES  ,12,/, 

OF  CONTROL  SYSTEM  INPUTS, 12,/, 

OF  RANDOM  NOISE  SOURCES  ,12,/, 

OF  DISPLAYED  OUTPUTS  ,12,/, 

29H  NO.  OF  DETERMINISTIC  INPUTS  ,12,/, 

29H  NO.  OF  TERMINAL  CONDITIONS  ,12,/, 

1H  ) 


29H  NO 
29H  NO 
29H  NO 


SPECIFY  DEL,  TO,  TEND,  AND  TEXT 
READ  (KIN,1200)  DEL,  TO,  TEND,  TEXT 
FORMAT  (8E10.0) 

TEND=IFIX((TEND-T0+0.0001 )/DEL)»DEL+T0 
WRITE  (KPTR,1225)  DEL,  TO,  TEND,  TEXT 
FORMAT  (25H  INTEGRATION  TIME  STEP  = ,F10.3,/, 
18H  INITIAL  TIME  = ,F10.3,/, 

18H  TERMINAL  TIME  = ,F10.3,/, 

18H  TIME  EXTENSION  = ,F10.3,/, 

1H  ) 

IDENTIFY  VARIABLES  FOR  OUTPUT 
CALL  OUTLET 
NSP=1 


DO  59  1=1,60 
IF  (LPRNT(I)  .GT. 
XMIN(I)=1.0E10 
XMAX(I)=-1 .0E10 
CONTINUE 


0)  NSP=NSPt-2 


SPECIFY  INTERNAL  TIME  BREAKS 
CALL  ipLET 
WRITE  (KPTR, 1310) 

FORMAT  (22H  INTERNAL  TIME  BREAKS  , 17H  INDEX 


CODE 


NDT) 
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1320 

65 

C 


C 

c 

70 

C 

C 

100 


150 

C 

205 


208 

C 

C 

C 

C 

210 

C 


215 

216 


C 

230 


DO  65  I=1,NC0DES 
IF  (NSTEP(I)  .EQ.  0) 
WRITE  (KPTR.1320)  I, 
FORMAT  (26X,  12,  3X, 
CONTINUE 

INITIALIZE  SOME  MORE 

QXT(1)=-1 .0 

QYT(1)=-1.0 

TCRO,  1)=9999.0 

NXP1=NX+1 

N=NX-fNU 

TIME=T0 

PTL=TIME 

XLsTIME 

TNEXTrTIME 

ND2=0 


GO  TO  65 

ICODES(I),  NSTEP(I) 
A4,  2X,  12) 


QUANTITIES 


MAIN  COMPUTATIONAL  LOOP  BEGINS  HERE 

START  BY  HANDLING  INTERNAL  AND  EXTERNAL  BREAKS 

CALL  UPDATE (TNEXT) 

IF  EITHER  THE  A OR  B MATRIX  HAS  CHANGED,  DISCRETIZE  THEM 
THEN  COMPUTE  NEW  FEEDBACK  GAINS 
IF  (IFLAG(1)+IFLAG(2)  .EQ.  0)  GO  TO  150 
CALL  DSCRT(N,AA.DEL,A,W0RK,4) 

CALL  EQUATE(BD.W0RK(1,NXP1),N,NU) 

CALL  EQUATE(SAV,A(1,NXP1),N,NU) 

CALL  GPFBGN(WORK) 

COMPUTE  EXP(A*T)  = (EXP(A»DEL))»*NPRED 
IF  (IFUG(1)+IFLAG(2)+IFLAG(10)  .EQ.  0)  GO  TO  210 
CALL  IDENT(N,EA,1.0) 

IF  (NPRED.EQ.O)  GO  TO  210 
CALL  EQUATE(EA,A,N,N) 

IF  (NPRED.EQ.l)  GO  TO  210 

DO  208  I=2,NPRED 

CALL  MMUL(A,EA,N,N.N,C0M2) 

CALL  EQUA TEC EA, COM2, N,N) 

CONTINUE 

DO  INCREMENTS  TO  MEANS 

P DISTRIBUTES  INCREMENT  BETWEEN  EST.  AND  ERR. 

SO  THAT  XBAR*"XHN“ENF 

S!sO  INCREMENT  SIGMA  MATRIX  AND  YBAR=DISPLAY  MEANS 
IF  (IFLAG(22) .EQ.O)  GO  TO  230 
P = 1 .0 
P=0 .0 

CALL  VMAT1(EA,XMINC,N,N,YBAR) 

DO  216  1=1, N 

XHN(I)=XHN  I)+YBAR(I)»P  ^ ^ 

ENF(I)=ENF(I)+XMINC(I)*(1 .0-P) 

XflAR(I)=XBAfl(I)+YBAR(I) 

SIGMA(Iij)=SIGMA(I,J)+XMINC(I)«XMINC(J)*(1 .0-P)**2 
CONTINUE 

CALL  VMAT1(C,XBAR,NY,N, YBAR) 

DO  INCREMENTS  TO  COVARIANCES 
IF  (IFLAG(23) .EQ.O)  GO  TO  250 
P=0.62 

DO  240  1=1, N , 

XCOVa,I)  = XCOV(I,IUXSINC  I)*»2 

ecovu.iUecov(i,iUxsinc(i  ••2*0 .0-P) 

SIGMA( I , I)  = SIGMA{ I , I)+XSINC( I )*»2*(  1 . 0-P) 
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DO  235  J 
YCOVCI.J 
YCOVa.I 
CONTINUE 


J=1.N 

ri,J)=6.i 

a,i)=P»: 


P»XSINC(I)*»2 


CALL  MMUL(EA,YC0V,N,N,N.C0M1) 

CALL  MAT6S(N,N,C0M1,EA,XHC0V) 

CALL  MAT3(NY,N,C,XC0V, YCOV) 

COMPUTE  MOTOR  NOISE 
DO  280  1=1, NU 
J=I+NX 

CHANGE  NEXT  STATEMENT  DEPENDING  ON  THE  CONTROLLER 
VUO(I)=VU(I)+PU(I)»ABS(XCOV(J, J)) 

W0(I+NW)=VU0(I) 

CONTINUE 

UPDATE  YNOM  FOR  NEW  XNOM 

IF  (IFLAG(25).EQ.1)  CALL  VMATKC, XNOM, NY, N, YNOM) 


COMPUTE  OBSERVATION  NOISE 
DO  235  1=1, NY 
YM=mR(I)-YNOM(I) 

YS=SQRT(YCOV(I.D) 

GTH(I)=XGAIN(TH(I).YM,YS)»»2 

VY0(I)=VY{I)+PY(I)*(YM»*2+YS«*2) 

VY0(I)=VY0(I)/GTH(I)/ATTN(I) 

CONTINUE 

PROPAGATE  COVARIANCES  ONE  TIME  STEP  AND  INCREMENT  TIME 
aLL  COVAR(WORK) 

TIME=TIME+DEL 

ADJUST  OUTPUTS  FOR  NOMINAL  PATH 
PROBLEM  DEPENDENT 
CALL  ADJNOM 

DO  OUTPUT  FOR  THIS  TIME  STEP 
CALL  INFORM(PTL,ND2,COM1,NSP) 

IF  THE  TIME  IS  NOT  EXPIRED  DO  ANOTHER  ITERATION 
IF  (TIME+0.0001  .LT.  TEND+TEXT)  GO  TO  70 


TIME  IS  EXPIRED.  DO 

CONTINUE 

XH=TEND+TEXT 

CALL  PRINTR(WORK,NSP) 

GO  TO  1 


DO  OUTPUT  AND  START  AGAIN 


RETURN  IF  AN  EOF  WAS  READ 
RETURN 


OOOOOo 
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nVCMP  - COMPUTATIONS  FOR  TIVAR 
INCLUDES 

1 - SUBROUTINE  GPFBGN 

2 - SUBROUTINE  TERMG 

- SUBROUTINE  COVAR 

- SUBROUTINE  FINAL 


i 


oKAL^PURPOsi°FlEDB»Clt  GAIN  COMPUTATION 
DIMENSION  X(1),  QXU(60) 

1 /COMMUN/  NCD,  PINT,  NSP,  ICODES(32)i  IFLAG(32) 

2 /INOU/  KIN,  KOUT,  KPTR  , , 

5 /MAIN1/  NDIM,  NDiAi,  COMKD 

■ /MAIN2/  Z(1)  ^ 

/C0MP5/  ADCL(I)  ^ ^ 

/INPTS/  BD(60),  DyM(60),  AD(1) 

/INPTX/  N,  NX,  rai,  NU,  AC(1) 

/INPTY/  NY,  Ch) 

/INPTW/  NW,  E(l) 

/WEIGHT/  QX^36)?^QY(30),  QR(30),  PSS(1) 
/GAINBK/  CGN(1) 


1 NM1=N-1 
NN=N»NDIM 
NNU=NU»NDIM 
NNW=NW«NDIM 

CALL^EQUA TE ( AD ( NNX ) , DUM , N , NU) 

C SKIP  STEADY  STATE  GAIN  COMPUTATION  UNLESS  NEWG  FLAG  IS  SET 

F (IFUG(21)  .EQ.  0)  GO  TO  100 

C STEADY  STATE  GAIN  COMPUTATION 

2 DO  5 1=1. NY 
C1=QY(I)*DEL/2.0 
DO  5 J=I,NN,NDIM 
X(J)=C1»6(J) 

5 CONTINUE 

CALL  MAT2A(NY,N,C,X,X) 

11=1 

X(II)=X(tl)+QX(I)*DEL/2.0 

II=II+NDIM1 

6 CONTINUE 

CALL  MAT3A(N.N,AD,X,Z)  , 

CALL  mscale(6xu,x(nnxKn,nu,del) 

11=1 

DO  9 I=1,NU  . , ^ 

Sll'^VSCALE(C0M1(II)  ,QXU(II)  ,N,-1 .0/QR(I+NU)) 

II=II+NDIM 

9 CONTINUE 

DO  10  1=1, N 

DO  10  J=I.NN,NDIM 

ADCL(J)=AD(J)^ 

Z(J)=Z(J)+X(J) 

X(J)=AD(J) 

10  CONTINUE 
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15 

C 


20 

1030 

1032 

1033 

C 

100 

C 

1 50 


1152 

170 

C 

200 


CALL  MAT5S(BD,C0M1,N,NU,N,p 
CALL  MAT6S(N,NU,QXU,C0M1,Z) 

CALL  TRANS1(N,X,X)  ^ ^ 

CALL  KFLTR(N,NU,X,BD,Z,QR(NU+1 ) ,PSS, DUM.CGN) 

CALL  MMUL(X,6uM,N,N,NU,Z) 

DO  15  1=1, NU 
n=fI-1)*NDIM  ^ 

C1=1.0/QR(NU+I) 

DO  15  J=I,NN,NDIM 
11=11+1 

CGN(J)=Z(II)+QXU(II)»C1 
COMUJ)  = -CGN(J) 

CONTINUE  , , 

CALL  MMULS(BD,C0M1,N,NU,N,ADCL) 

GET  EQUIVALENT  CONTINUOUS  GAINS 

C1=-DEL 

DO  20  L=1,2 

CALL  EQUATE(X,AC,N.N) 

CALL  EQUATEU(nX+1  ).COM1,NU.N) 

CALL  DSCRT(N,X,DEL,Z,COM1,J^) 

CALL  GMINV(N,N,C0M1,Z,MR, 1) 

CALL  MMUL(CGN.Z,NU,N,N,C0M1) 

IF  (L.EQ.2)  CUDEL 

CALL  MSCALE(C0M1,C0M1,NU,N,C1) 

CONTINUE 

FORMAT^ (U3H’ COMPUTED  STEADY  STATE  DISCRETE  CTRL  GAINS  ), 
CALL  MATPRT(NU,N.CGN,5HDGAIN) 

F0RMAT^?35H’ EQUIVALENT  CONTINUOUS  GAINS  LX,LU  ) 

CALL  MATPRT(NU.N,COM1,5HCGAIN) 

CALL  GM1NV(NU,NU,C0MUNNX)  ,Z(NNX),MR,  1) 

CALL  MMUL ( Z( NNX ) , COM  1 , NU , NU , NX , Z ) 

FORMAT^ (31H’ EQUIVALENT  CONTINUOUS  LOPT,TN  ) 

CALL  MATPRT(NU,N,Z,5HL*,TN) 

IF  GAINS  WERE  NOT  READ  IN,  COMPUTE  TERMINAL  GAINS 
IF  (IFLAG(19)  .EQ.  1)  GO  TO  150 
IF  (iFLAGns)  .EQ.  1)  GO  TO  170 
CALL  TERMG(X,QXU) 

GO  TO  200 

DGAINS  WERE  READ  IN 
CALL  EQUATE(C0M1,AC,NX,N) 

CALL  MSCALE(COMUNX+1 ) ,CGN,NU.N,-1 .0) 

CALL  DSCRT(N,COM1,DEL,CGN,Z,^) 

CALL  MMUL ( COM KNX+1  ) , Z,  NU,  N,  N,  CGN) 

CALL  MSCALE(CGN.CGN,NU,N,-1 .O/DEL) 

FORMAT^ (37h’ EQUIVALENT  DISCRETE  GAINS  GENERATED  ) 

CALL  MATPRT(NU,N,CGN,5HDGAIN) 

CONTINUE 
KILL  TCOST  FLAG 

CALL  MMUL(BD,CGN(NNX) ,N,NU,NU,E(NNW+1 )) 

DO  205  I=NNX,NN,NDIM 

II=I+NM1 

DO  205  J=I,II 

K=J-NX*NDIM 

L=K+NNW 

DUM(K)=AD(J)  , , 

AD(J)=AD(J)-E(L) 


II 


I 


44 
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E(L)=E(L)/DEL 
205  CONTINUE 
RETURN 

END 


'J’oo 
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r 

i 

t 


c 


1 


1 

2 

5 

8 

9 

A 

B 

E 

F 

G 

J 

1 

C 

C 


8 


10 


15 


16 


17 

C 

C1800 


SUBROUTINE  TERMG(X, QXU) 

COMPUTES  CONTROL  GAINS  BASED  ON  TERMINAL  COSTS  (FORWARD  TIME) 
DIMENSION 

X(1),  QXUd),  BR(60),  FM(225),  PHN(225) 


COMMON 

/COMMON/ 

/INOU/ 

/MAIN1/ 

/MAIN2/ 

/COMP5/ 

/INPTS/ 

/INPTX/ 

/INPTY/ 

/TIMES/ 

/WEIGHT/ 

/TCOST/ 

/GAINS K/ 


NCD,  PINT,  NSP,  ICODESC32),  IFLAG(32) 
KIN,  KOUT,  KPTR 
NDIM,  NDIM1,  COMI(l) 

C0M2(1) 


ADCL(1) 
BD(60), 


AD(1) 


SAV(60) 

N,  NX.NX1,  NU 
Nir,  ch) 

TIME.  DEL.  TO,  TEND,  TEXT,  TD,  NPRED 
QX(30),  QY(30),  QR(30),  PSS(1) 
QXT(30),  QYT(30),  PINV(1) 

CGN(1) 


NN=N*NDIM 

RETURN  IF  NO  TERMINAL  COSTS 

ff  (QXT(1)+QYT(1)  .LE.  -1.5)  RETURN 

CHECK  TIME  TO  GO 

NGO= ( TEND-TIME+0 . 000 1 ) /DEL-1 

IF  (NGO)  150,50,5 

DO  THIS  PART  IF  NGO  IS  POSITIVE 

AND  IF  QXT  OR  QYT  HAVE  CHANGED 

IF  aFLAG(27)+IFLAG(28)  .EQ.  0)  GO  TO  50 

IF  (TIME  .GT.  TO+0.0001)  GO  TO  8 

NGO=NGO-NPRED 

CALL  EQUATE (X,ADCL,N,N) 

CALL  MMUL(PSS,BD,N,N,NU,C0M1) 

CALL  MAT2A(N,NU,BD,C0M1,C0Mt) 


11=1 

DO  10  1=1, NU 

C0M1(II)=C0M1(II)+QR(I)»DEL 

II=II+NDIM1 

CONTINUE 

CALL  GMINV(NU,NU,COM1,COM2,MR, 1) 
CALL  MMULUd,COM2,N,n0,NU,BR) 
CALL  MAT6(N,NU,BR,BD,PHN) 

CALL  DINTEG(N,X,PHN,PINV,NG0) 
TPT=TIME 


CIrO.O 

DO  15  1=1, NY 

ff  (QYT(I)  .GT.  -0.5)  C1=QYT(I) 

DO  15  J=I,NN,NDIM 
X(J)=C1»C(J) 

CONTINUE 

CALL  MAT2A(NY,N,C,X,X) 

11=1 

DO  17  1=1, N 

DO  16  J=I,NN,NDIM 

X(J)=X(J)-PSS(J) 

CONTINUE 

JF  (QXT(I)  .GT.  -0.5)  X(II)=X(II)+QXT(I) 

II=II+NDIM1 

CONTINUE 

WRITE  (KOUT, 1800)  TEND 

FORMAT  (27H  TERMINAL  CONDITION  AT  TIME,F10.3) 
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C CALL  MATI0(X,N,N,3) 

CALL  GMINV(N,N,X,FM,MR, 1) 

DO  25  1=1, N 

DO  25  J=I,NN,NDIM  , ^ 

00M2(J)=FM(J)+PINV(J) 

25  CONTINUE  „ . 

CALL  GMINV(N,N,C0M2,C0M1,MR,  1) 

CALL  MMUL(C0M1,PHN,N.N.N,C0M2) 

CALL  MAT2A(N.N,PHN.c6M2,C0Mn 
CALL  MSCALE(BR,BR,N,NU,-1 .0) 

GO  TO  100 

C UPDATE  PINVERSE 

C DO  THIS  PART  IF  NGO  IS  ZERO 

C OR  IF  NEITHER  QXT  OR  QYT  HAVE  CHANGED 

50  IF  (TIME-TSS)  62,62,60  ^ 

60  CALL  MMUL(ADCL,FM,N,N,N,C0M1) 

CALL  MAT2(N,NU,BR,BD.FM) 

CALL  MAT6S(N,N,C0M1,ADCL,FM) 

62  DO  70  1=1, N 

DO  70  J=I,NN,NDIM  ^ 

C0M2(J)=FM(J)+PINV(J) 

calPgminv(n,n,com2,comi,mr,  1) 

CALL  MMUL( COM 1,PHN,N.N.N, COM2) 

CALL  MAT2A(N,N,PHN,C0M2,C0M1) 

C GET  NEW  CONTROL  GAINS 

C DO  THIS  PART  UNLESS  NGO  IS  NEGATIVE 

100  DO  110  1=1, N 

DO  110  J=I,NN.NDIM^ 

X(J)=PSS(J;+C0M1(J) 

110  CONTINUE 

C WRITE  (KOUT.1900)  TIME 

C1900  FORMAT  (26H  DIFFERENCE  MATRIX  AT  TIME.F10.3) 

C CALL  MATI0(C0M1,N,N,3) 

CALL  MMUL(X,BD,N,N,NU,COM1) 

CALL  MAT2A(N,NU,BD, COM 1, COM2) 

11=1 

DO  112  I = 1.NU^  ^ 

COM2(II)=c6m2(II)+QR(I)*DEL 
n=II+NDIM1 
112  CONTINUE 

CALL  GMINV(NU,NU,C0M2,X,MR, 1) 

CALL  TRANS2(N.NU,QXU,C0M2) 

CALL  MAT5AS( COM  1 , AD , NU, N, N, COM2 ) 

CALL  MMUL(X,C0M2,NU,NU,N,CGN) 

IF  (TIME.LT.TPT)  GO  TO  150 
115  TPT=TIME+PINT 

WRITE  (KOUT.1950)  TIME 

1950  FORMAT  (21H  CONTROL  GAIN  AT  TIME.F10.3) 

CALL  MATIO(CGN,NU,N,3) 

150  RETURN 


END 
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SUBROUTINE  C0VAR(X5) 

PROPAGATES  COVARIANCES  ONE  TIME  STEP 


DIMENSION 
X5(l),  DELF(5,30) 


COMMON 

1 /COMMUN/ 

2 /INOU/ 

4 /MAIN1/ 

5 /MAIN2/ 

9 /INPTS/ 

A /INPTX/ 

B /INPTY/ 

C /INPTW/ 

D /INPTWD/ 
E /TIMES/ 

H /NOISE/ 

J /GAINBK/ 
L /COVX/ 

M /COVY/ 

N /COVXHE/ 
0 /COVXH/ 

P /COVEF/ 

Q /NOMNL/ 


NCODES,  PINTVL,  NSP,  ICODES(32),  IFLAG(32) 
KIN,  KOUT,  KPTR,  KPUNCH,  KDISK 
NDIM,  NDIM1,  X1(1) 

X3(i5 

BD(60),  SAV(60),  A(1 ) 

N,  NX,  NX1,  NU 
Ni,  Cl) 

NW,  EU) 

NWD,  F(1) 

TIME,  DEL,  TO,  TEND,  TEXT,  TD,  NPRED,  EA( 1 ) 

vu(36),  vy(3o),  wo(3o),  sigma(i) 

CGN( 1 ) 

XBARC30),  XCOV(1) 

YBAR(30),  YCOV(I) 

DUM(30),  XHECOVd) 

XHN(30),  XHCOVd) 

ENF(30),  ENP(30),  ECOVd) 

NTF,  XNOMd) 


IF  (TIME  .GT.  T0+0.001)  GO  TO  8 

INITIALIZATION 
NWU=NW-»-NU 
hM1=N-1 
NN=N«NDIM 
DO  1 K=1,NWD 

DELF(K,1)=FDOT(K,TIME)«DEL 

CONTINUE 

TCOR=1 ,0 

TCOR=SQRT(TCOR/DEL) 

IF  (NPRED. EQ.O)  GO  TO  8 

INITIALIZE  FOR  THE  TIME  DELAY 
DO  2 Izl, NPRED 
CALL  FINAL(X5) 

IFLAG(26)=0 
TIME=TIME+DEL 
DO  2 K=1,NWD 

DELF(K,I+1 )=FDOT(K,TIME)»DEL 
CONTINUE 

INTEGRATE  THE  DRIVING  AND  MOTOR  NOISES 
11=1 

DO  9 1=1, NWU 

™p=wo(i)«del 

CALL  VSCALE(X5(II) ,E(II) ,N,TMP) 

II=II+NDIM 

CONTINUE 

CALL  MAT2(N,NWU,E,X5,X5) 

COMPUTE  COVARIANCES  AND  ESTIMATOR  GAIN 
DO  10  1=1, NY 
TMP= 10000.0 

IF  (VY(I)  .GE.  1.0E-5)  TMP=DEL/VY(I) 

DO  10  J=I,NN,NDIM 

X3(J)=C(J)*TMP 

CONTINUE 
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CALL  MAT2A(NY,N,C,X3,X3) 

CALL  MMUL(S1GMA,X3,N,N,N,XC0V) 

DO  11  I=1,NN,NDIM1 
XC0V(I)rXC0V(I)+1 .0 
CONTINUE 

CALL  GMINV(N,N,XCOV ,iCGV,MR,  1) 

CALL  MMULCYCOV, SIGMA,  N.N.N.XCOV) 

CALL  MMUL(XC0V,X3,N,N,N,X1) 

CALL  D1AG2(N,YC0V,X1,-1 .0, 1.0) 

CALL  MMUL(EA,X1,N,N,N,X3j 
CALL  MAT2(N,N,XC0V,X1,X1) 

CALL  MMUL(A,X1,N,N,N,XCOV) 

CALL  MAT6S(N,N,XCOV,A,X5) 

CALL  MMULU,YC0V,N,N,N,X1) 

CALL  MMULCEA, SIGMA, N,N,N,YC0V) 

CALL  MMUL(X1,EC0V,N,N,N,XC0V) 

CALL  MAT6S(N,N,XC0V,X1,X5) 

CALL  VMAT2(XHN,X3,ENF,N,N,pUM) 

CALL  VMAT1(X1,ENF,N,N,XBAR) 

DO  15  1=1, N 
ENF(I)=XBAR(I) 

DO  15  J=I,NN,NDIM 

SIGMA { J ) = ECOV ( J ) -S IGMA ( J ) 

EC0V(J)=X5(J) 

CONTINUE 

CALL  MMUL(X3, SIGMA, N,N,N,X5) 

NNX=NX*NDIM 
DO  20  1=1, N 
IMUI-1 

DO  20  Jrl.NN.NDIM 

XCOV(J)=XHECOV(J)+(YCOV(J)+X5(J))/2.0 

YCOy(J)=A(J) 

IF  (J  .GT.  NNX)  GO  TO  20 

YCOV(J)  = YCOV(J)-DOT3(NU,BD(I)  ,CGN(J-IMD) 

XHEC0V(J)=XHEC0V(J)+X5(J) 

CALL  MMUL(YCOV,XHECOV,N,N,N,X5) 

CALL  MAT5(X5,X1,N,N,N,XHECOV) 

CALL  MAT5(Xc6v,X3,N,N,N,X5) 

CALL  MMUL(X1, SIGMA, N,N,N,X3) 

CALL  MAT2(N,N.X3, XI, SIGMA) 

CALL  MMUL(EA,ECOV,N,N,N,X3) 

CALL  VMAT1(YCOV,DUM,N,N,XHN) 

DO  30  II=1,NN,NDIM1 
K=(II+NDIM)/NDIM1 
DUM(K)=0.0 
I = II 

DO  30  J = II,NN,NDIM 

XHCOV{J)=XHCOV(J)+X5(J)+X5(I) 

SIGMA(J)=EC0V(J)-SIGMA(J) 

X1(J)=XHEC0V(J)+X3(J)/2.0 

X1(I)=XHECOV(I)+X3(I)/2.0 

SIGMA(I)=SIGMA(J) 

XHCOV(l)=XHCOV(J) 

1 = 1+1 
CONTINUE 

CALL  MMUL(YCOV,XHCOV,N,N,N,X3) 

CALL  MAT2(N,N,X3, YCOV,XHCOV) 

CALL  MAT5(X1,EA,N,N,N,X3) 

DO  W 11=1 ,NN,NDIM1 
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I=II 

DO  40  J=II,NN,NDIM 

XCOV(J)=XHCOV(J)+X3(J)+X3(I) 

XC0V(I)=XC0V(J) 

1=1+1 

CONTINUE 


n=i 

DO  50  1=1, NWD 
TMP=TC0R»DELF(I,1) 
CALL  VSCALE(X5(II), 
II=II+NDIM 
CONTINUE 


F(II),NX,TMP) 


CALL  VMAT2(ENF,F,DELF,NX,NWD,ENF) 

CALL  MAT6S(NX,NWp,X5,X5, SIGMA) 

IF  (NPRED  .EQ.  0)  GO  TO  70 
11=1 

DO  51  1=1, NWU 
TMP=SQRT(WO(I)«DEL) 

CALL  VSCALE(X1(II),E(II),N,TMP) 

II=II+NDIM 

CONTINUE 

DO  60  1=1, NPRED 

CALL  MAT6S(N,NWU,X1,X1,XCOV) 

IF  (I  .EQ.  NPRED)  GO  TO  53 
CALL  MMUL(A,X1,N,N,NWU,X3) 

CALL  EQUATE(X1,X3,N,NWU) 

DO  55  K=1,NWD 

delf(k,i)=delf(k,i+i ) 

CONTINUE 

F (I  .GT.  1)  CALL  VMAT1(A,ENP,NX,NX,DUM) 
CALL  VMAT2(DUM,F,DELF(1 ,1) ,NX,NWD,ENP) 
CONTINUE 

CALL  VMAT2(ENP,EA,ENF,N,N,DUM) 

CALL  FINAL(X5) 

DO  75  1=1, N 

XBAR(I )=XHN( I )+DUM( I )+XNOM( I ) 

CONTINUE 

CALL  MAT3(NY,N,C,XCOV,YCOV) 

CALL  VMAT1(C,XBAfi,NY,N,YBAR) 

DO  80  K = 1,NWD 

DELF(K,NPRED+1 )=FDOT(K,TIME+DEL)»DEL 
CONTINUE 

RETURN 

EMD 
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SUBROUTINE  FINAL(X) 

COMPUTES  NOMINAL  PATH  BASED  ON  TERMINAL  CONDITION 
CONSTRAINED  BY  MINIMIZING  CTRL  RATE 

DIMENSION 
1 Q0(30),  X(l) 


COMMON 

/COMMUN/  NCODES 
/MAIN1/  NDIM, 


/MAIN2/ 

/C0MP1/ 

/COMP3/ 

/INPTS/ 

/INPTX/ 

/INPTY/ 

/TIMES/ 


NCODES,  PINTVL,  NSP,  ICODE 
NDIM,  NDIM1,  COMKI) 

Y(1) 

BD(60),  SAV(60),  AD(  1 ) 

N,  NX,  NX1,  NU,  AA(1) 

NY.  C(1) 

TIME.  DEL,  TO,  TEND,  TEXT, 


ICODES(32),  IFLAG(32) 


/ xxricivj/  xju'icif  uau*  xw,  xxiniy,  x x*aT  , TD, 

/WEIGHT/  QX(30),  QY(30),  QR(3oj 
/NOMNL/  NT,  XN0M(30;,  YN0M(30),  P(l) 

P(1)=9999  IS  A FLAG  TO  AVOID  COMPUTATION 
IF  (P(l5  .EQ.  9999.0)  RETURN 
NN=N»NDIM 

NGO=(TEND-TIME+Q.0001)/DEL 

IF  (NGO  .LE.  NT)  GO  TO  50 

ASSUMES  SYSTEM  IS  OUTPUT  CONTROLLABLE 

IF  (IFLAG(26).EQ.O)  GO  TO  50 
A NEW  TCR  MATRIX  WAS  READ  IN 
NNX=NX»NDIM 
11=1 

DO  6 1=1, NU 
C1=1 .0/(OR(I)*DEL) 

L=II+N-1 
DO  5 J=II,L 
C0M1(J)=BD(J)»C1 
AD(J+NNX)=SAV(J) 

CONTINUE 

II=II+NDIM 

CONTINUE 

CALL  MAT2(N,NU,COM1,BD,Z) 

CALL  EQUATE(ADI,AD,N,N) 

CALL  DINTEG(N,ADI,Z,X,NGO) 

CALL  MAT3(NT,N,P,X,C0M1) 

CALL  GMINV(NT,NT,c6m1,X,MR, 1) 

CALL  MMUL(P,Z,NT,N,N,COM1) 

CALL  VMAT2(P(NN+1 ),C0M1,XN0M,NT,N,Q0) 
CALL  VMAT1(X,Q0,NT,NT, YNOM) 

CALL  TRANS2(NT,N,C0M1,Y) 

CALL  VMAT1(Y,YN0M,N,NT,Q0) 

CALL  TRANS1(N,AD,COM1) 

CALL  GMINV(N,N,c6m1,ADI,MR, 1) 

CALL  EQUATE(AD(NNX+1 ),SAV,N,NU) 

UPDATE  XNOM  AND  YNOM 

IF  (NGO  .GT.  0)  GO  TO  75 

DO  70  1=1,  N 

Q0(I)=0.0 

CONTINUE 

11=1 

CALL  VMAT1(ADI,Q0,N,N, YNOM) 

DO  77  1=1, NU 

C0M1(I)=-D0T(N,BD(II) ,YNOM)/(QR(I)»DEL) 
II=II+NDIM 


NPRED 
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77 


CONTINUE 


CALL  VMAT1(AD,XNOM,N,NX,QO) 
CALL  VMAT2(Q0,SAV,XN0M(NX+1 
CALL  VMAT2(OO.BD,iOM1,N,Ny, 
CALL  VSCALE(q6,YNOM,N,1.0) 
CALL  VMAT1(C,XN0M,NY,N,YN0M) 


N.NU.QO) 


RETURN 

END 
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nvio  - I/O  Fo;.  tivar 

INCLUDES 

1 - SUBROUTINE  INTLET 

2 - SUBROUTINE  UPDATE 

3 - SUBROUTINE  OUTLET 

4 - SUBROUTINE  INFOriM 

5 - SUBROUTINE  PRINTR 


SUBROUTINE  INTLET 

SPECIFY  INTERNAL  UPDATES 

NSTEP(I)  REFERS  TO  THE  ITH  TYPE  OF  UPDATE 

NSTEP(I)=0  NO  INTERNAL  UPDATE 

NSTEP(I)=NDT  NDT  DELS  BETWEEN  SUCCESSIVE  INTERNAL  UPDATES 


1 

1 

2 


COMMON 


/COMMUN/  NCODES,  PINTVL,  NSP,  IC0DESp2),  IFLAG(32), 
NSTEP(32),  LPRNT(60),  XMAX(60),  XMIN(60) 
/INOU/  KIN,  KOUT,  KPTR,  KPUNCH,  KDISK 


C READ  NSTEP(I)  FROM  SEVERAL  CARDS 

100  READ  (KIN, 1100)  (NSTEP(I),  1=1, NCODES) 

1100  FORMAT  (1615) 

RETURN 

END 
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c 


1 

1 

2 

4 

A 

B 

C 

D 

E 

F 

G 

H 

I 

J 

K 

Q 


1 


100 

C 

1 10 

c 

120 

1030 


C 

130 

135 


C 

140 


k 


c 

150 

C 


160 

C 


SUBROUTINE  UPDATE(TNEXT) 
PERFORMS  EXTERNAL  UPDATES 


COMMON 
/COMMON/ 

/INOU/ 

/MAIN1/ 

/INPTX/ 

/INPTY/ 

/INPTW/ 

/INPTWD/ 

/TIMES/ 

/WEIGHT/ 

/TCOST/ 

/NOISE/ 

/RATIOS/ 

/GAINBK/ 

/INC RE/ 

/NOMNL/ 

DATA 

LEND,  12,  LTIME,  PI  /3HEND,  1HZ,  4HTIME, 


NCODES,  PINTVL,  NSP,  IC0DES(32) 
N STEP (32) 

KIN,  KOUT,  KPTR,  KPUNCH,  KDISK 

NDIM,  NDIM1 

N,  NX,  NX1,  NU,  AA(1) 

NY,  Ch) 

NW,  E(1) 

NWD,  F(1) 

TIME.  DEL,  TO,  TEND.  TEXT,  TD, 

qx(36),  qy(3oJ.  qr(3o) 

QXT(30),  QYT(30 

vuoho),  vYoho).  wo(30)  , ^ 

PU(30),  vu(30),  fY(30),  VY(30), 

GAINS( 1 ) 

XMINC(30),  XSINC(30),  ATTN(30) 
NTF,  XNOM(30),  YNOM(30),  TCR(1) 


, IFLAGC32), 


NPRED 

TH(30),  GTH(30) 


3.14159/ 


NN=N+1 

NNX=NX+1 

NNN=NX»NDIM+1 

TMAX=TEND+TEXT+0 . 1 

TDMAX=29*DEL 

TAKE  CARE  OF  THE  INTERNAL  UPDATES 
DO  100  1=1, NCODES 
IFUG(I)  = 0 

IF  (NSTEP(I) .EQ.O)  GO  TO  100 
TTME=TIME-TO 

IF  (M0D(INT(TTME/DEL+0.1),NSTEP(D)  .EQ.  0)  CALL  INTNEW(I) 
CONTINUE 

RETURN  IF  NOT  TIME  FOR  THE  NEXT  EXTERNAL  UPDATE 
IF  (TIME+0.001  .LT.  TNEXT)  GO  TO  220 

SPECIFY  THE  TYPE  OF  EXTERNAL  UPDATE  (AND  TIME  IF  LTIME) 
READ  (KIN, 1030)  IDEN,  BRKT 
FORMAT  (A5,5X,E10.0) 

IF  (EOF(KIN))  135,130 

CHECK  FOR  LEND 

IF  ( IDEN. NE. LEND)  GO  TO  140 

TNEXTsTMAX 

GO  TO  110 

CHECK  FOR  LTIME 

IF  (IDEN. NE. LTIME)  GO  TO  150 

TNEXTrBRKT 

GO  TO  110 

CHECK  FOR  LZ  (NULL  UPDATE)  - FOR  HRA  COMPATABILITY 
IF  (IDEN.EQ.LZ)  GO  TO  120 

SEARCH  THROUGH  THE  UPDATE  CODES,  ICODE(KEY) 

DO  160  KEY=1, NCODES 

IF  (IDEN.EQ.ICODES(KEY))  GO  TO  170 

CONTINUE 

CODE  WAS  ILLEGAL 
WRITE  (KOUT, 1165)  IDEN 
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1165 

C 

170 

1175 

C 

1 

2 

3 

4 

5 

C 

6 

7 

8 

9 

C 

10 


C 

11 

12 

13 
200 

14 

15 
201 

C 

16 

17 

C 

18 


FORMAT  (1H  ,5HC0DE  ,A5,HH  IS  ILLEGAL) 
CALL  EHT 


DO  THE  SPECIFIED  EXTERNAL  UPDATE 
IFUG(KEY)=1 
10=2 


FORMAT^f A2§H^Ix‘^ERNAL’UPM^  AT  TIME  .F8.3,4X,4HCODE.2X.A5) 

SYSTEM  DYNAMICS  - A,  B,  C,  D,  E 
CALL  MATI0(AA,NX,NX,I0) 

GO  TO  120 

CALL  MATI0(AA(NNN),NX,NU,I0) 

GO  TO  120 

CALL  MATI0(C,NY,NX,I0) 

GO  TO  120 

CALL  MATIO(C(NNN) ,NY,NU,I0) 

GO  TO  120 

CALL  MATIO(E,NX,NW,IO) 

GO  TO  120 


COST  FUNCTIONAL  - QX,  QY,  QU,  QR 
CALL  VECTIO(QX,NX,IO) 

GO  TO  120 

CALL  VECTIO(QY,NY,IO) 

GO  TO  120 

CALL  VECTIO(QX(NNX),NU,IO) 

GO  TO  120 

CALL  VECTIO(QR,NU,IO) 

GO  TO  120 

TIME  DEUY  - TD 
CALL  VECTIO(TD,1,1) 

NPRED=IFIX((TD+0.0001)/DEL) 

TD=DEL»NPRED 
CALL  VECTIO(TD, 1,3) 

GO  TO  120 

NOISES  - PN,  MNA,  MNR,  SNA,  SNR 
CALL  VECTIONo,NW,IO) 

GO  TO  120 

CALL  VECTI0(VU,NU,I0) 

GO  TO  120 

CALL  VECTI0(PU,NU,I0) 

DO  200  1=1, NU 

PU(I)=PI*l6.0«*(PU(I)/10.0) 

GO  TO  120 

CALL  VECTIO(VY,NY,IO) 

GO  TO  120 

CALL  VECTIO(PY,NY,IO) 

DO  201  1=1. NY 

PY(I)=PI»i6.0«»(PY(I)/10.0) 

GO  TO  120 

THRESHOLDS  AND  ATTENTION  - TH,  ATTN 
CALL  VECTIO(TH,NY,IO) 

GO  TO  120 

CALL  VECTIO(ATTN,NY,IO) 

GO  TO  120 

CONTINUOUS  AND  DISCRETE  CTRL  GAINS  - CGAIN,  DGAIN 
CALL  MATIO(GAINS,NU,N,IO) 

GO  TO  120 
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19  CALL  MATIO(GAINS,NU,N,IO) 

GO  TO  120 

C PRINT  INTERVAL  - PNTVL 

20  CALL  VECTI0(PNTVL,1,I0) 

GO  TO  120 

C COMPUTE  NEW  GAINS  - NEWG 

21  GO  TO  120 

C INCREMENTS  TO  STATES  - XMINC,  XSINC 

22  CALL  VECTIO(XMINC,N,IO) 

GO  TO  120 

23  CALL  VECTIO(XSINC,N,IO) 

GO  TO  120 

C DETERMINISTIC  INPUT  (F  MATRIX)  - F 

2M  CALL  MATIO(F,NX,NWD,IO) 

GO  TO  120 

C NOMINAL  PATH  VECTOR  - XNOM 

25  CALL  VECTIO(XNOM,N,IO) 

GO  TO  120 

C TERMINAL  CONDITION  MATRIX  - TCR 

26  CALL  MATIO(TCR,NTF,NN,IO) 

GO  TO  120 

C CALL  AN  INTERNAL  UPDATE  - INT 

27  READ  (KIN. 1210)  KEY 

1210  FORMAT  (12) 

WRITE  (K0UT.1220)  KEY 
1220  FORMAT  (17H  INTNEW  PART  NO.  ,13) 

CALL  INTNEW(KEY) 

GO  TO  120 

C TERMINAL  COSTS  - QXT,  QYT 

28  CALL  VECTIO(QXT,NX, 10) 

GO  TO  120 

29  CALL  VECTIO(QYT,NY,IO) 

GO  TO  120 

C DUMMY  UPDATE  - DUMMY 

30  CONTINUE 
GO  TO  120 

C NO  MORE  EXTERNAL  UPDATES  AT  THIS  TIME 

220  CONTINUE 

C GAIN  COMPUTATION  REQUEST  IS  CANCELLED  BY  READ-INS 

IF  (IFLAG(18).EQ.1  .OR.  IFLAG( 19) .EQ. 1 ) IFLAG(21)=0 
RETURN 

END 


I 
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SUBROUTINE  OUTLET 

fpPNTM  2of STAtIs-  LPRNT(21-1»0)  OUTPUTS:  LPRNT(41-60)  CONTROLS 

LPRlSUi^l  PRINTfllPRNTU  PRINT  & PLOi;  LPRNT(I)-3  PLOT 

1 /COMMUN/  NCODES,  PINTVL,  NSP.  ICODES(32),  IFLAG(32), 

1 NSTEP(32),  LPRNT(60) 

fTTkl  VATIT  IfDTR 


2 /INOU/  KIN,  KOUT,  KPTR 
4 /MAIN1/  NDIM 


C 

100 

1100 


READ  LPRNTd 
READ  (KIN.n 
FORMAT  (6011 
RETURN 

END 


) FROM  ONE  CARD 

00)  Tlprnt(i), 


1=1,60) 
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SUBROUTINE  INF0RM( PTL, NPTL , BUF , KDIM) 

DO  OUTPUT  FOR  A SINGLE  TIME  STEP 

DIMENSION  BUF(KDIM) 

COMMON 

/COMMUN/  NCODES,  PINTVL,  NSP,  IC0DES(32).  IFLAG(32), 

NSTEP(32),  LPRNTX(26),  LPRNTY(20),  LPRNTU(20), 
XMAX(20),  YMAX(20),  UMAX(20), 

XMIN(20),  YMIN(20),  UMIN(20) 

/INOU/  KIN,  KOUT,  KPTR,  KPUNCH,  KDISK 
/MAIN1/  NDIM,  NDIM1 

N.  NX 

TIME 


/INPTX/ 

/TIMES/ 

/COVX/ 

/COVY/ 


XBAR(30),  XCOVd) 
YBAR(30),  YCOVd) 


LL=1 

CHECK  IF  TIME  FOR  A PRINTOUT 
IF  (TIME-PTL+0.001  .LT.  PINTVL)  RETURN 

DO  A PRINTOUT 
PTL=TIME 
NPTL=NPTL+1 
BUF(1)=TIME 

DO  STATES 
DO  10  L=2,4.2 
DO  10  1=1,26 
LPIrLPRNTXd) 

IF  (LPI.NE.L  .AND.  LPI.NE.L-1)  GO  TO  10 
LULL+2 

BUF(LL-1 )=XBAR(I) 

II=I*NDIM1-NDIM 

BUF(LL)=SQRT(ABS(XCOV(II))) 


IF  (LPI.EQ.I)  GO 
C1=XBAR(iy+BUF(LL 
C2=XBAR(I)-BUF(LL, 
IF  (CI.GE.XMAXdl, 
IF  (C2.LE.XMIN(I). 
CONTINUE 


ro  10 


XMAX(I)=C1 

XMIN(I)=C2 


DO  OUTPUTS 
DO  20  L=2,4.2 
DO  20  1=1,26 
LPI=LPRNTY(I) 

IF  (LPI.NE.L  .AND.  LPI.NE.L-1)  GO  TO  20 
LL=LL+2 

BUFCLL-1 )=YBAR(I) 

II=I*NDIM1-NDIM 

BUF(LL)=SQRT(ABS(YCOV(II))) 

IF  (LPI.EQ.I)  GO  TO  20 
C1=YBAR(I)+BUF(LL) 

C2=YBAR(I)-BUF(LL 


IF  (CI.GE.YMAX(I), 
IF  (C2.LE.YMINa), 
CONTINUE 


YMAX(I): 

YMIN(I): 


Cl 

C2 


DO  CONTROLS 
DO  30  L=2,4,2 
DO  30  1=1,26 
LPI=LPRNTU(I) 

IF  (LPI.NE.L  .AND.  LPI.NE.L-1)  GO  TO  30 

J=I+NX 

LL=LL+2 

BUF(LL-1 )=XBAR(J) 
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30 


JJ=J«NDIM1-NDIM 
BUF(LL)=SQRT( ABS(XC0V( JJ) ) ) 
IF  aPI.EQ.D  GO  TO  30 


C1=XBAR(J)+ 

C2=XBAR(J)- 


C2=XBAR 
IF  (Cl. 


+BUF(LL) 

-BUF(LL) 


C1.GE.UMAX( 

C2.LE.UMIN( 


CONTINUE 


C WRITE  THE  BUFFER  ON  THE  DISK  AND  RETURN 

WRITE  (KDISK)  BUF 
RETURN 
END 


I 


I 
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C 


9 

10 


30 

C 

35 

36 


45 

40 

C 

50 


1 


1 

1 

2 

3 

I 

A 


1 

1 

2 

2 


SUBROUTINE  PRINTR( BUF , KDIM) 

DO  ALL  THE  OUTPUT  AT  THE  END  OF  A RUN 

DIMENSION 

BUFCKDIM),  TITLE(3),  LET(10,3) 

COMMON 

/COMMUN/  NCODES.  PINTVL,  NSP,  ICODES(32),  IFLAG(32), 
NSTEP(32),  LPROT(2o!3).  XMAX(60),  XMIN(60) 

/INOU/  KIN,  KOUT,  KPTR,  KPONCH,  KDISK 

/PLOT1/  NV,  NH,  NCPW,  LW,  XL,  XH,  YL,  YH,  NXES,  NDIR,  1ST, 
NGLV,  NGLH,  BSYM,  GSYM,  PSYM,  ND1 , ND2,  NOUT 
/MAIN1/  NDIM,  NDIMU  GRAPH(I) 

/INPTX/  STORE(l) 

DATA 

PS1,  PS2,  LMEAN,  LSD 

/1HM.  1HS,  4HMEAN,  8HSTD  DEV/, 

(TIT^Ed),  1=1,3) 

/8H  STATE  ,8H  OUTPUT  ,8HC0NTR0L  / 

IBEG=ND2+1 
REWIND  KDISK 

DO  10  1=1, ND2 
READ  (KDISK)  BUF 
II=I 

DO  9 L=1,KDIM 
STORE(II)=BUF(L) 

II=II+ND2 

CONTINUE 

DO  90  1=1,3 
M=0 


J=1,M) 


DO  30  L=1.20 
IQ=LPRNT(L.I) 

IF  (IQ.EQ.O  .OR.  IQ.EQ.3)  GO  TO  30 
M=M+1 

LET(M,I)=L 
CONTINUE 

IF  (M.EQ.O)  GO  TO  50 

PRINT  OPTIONS 
CALL  PAGEFD(KPTR.I) 

WRITE  (KPTR, 35)  (TITLE(I),  LET(J,I), 

FORMAT  (1H  ,9X.5(5X,A8,I4,5X)) 

WRITE  (KPTR, 36)  ( LMEAN, LSD, J=1 ,M) 

FORMAT  (1H  ,2X,4HTIME,2X,5(5X,A4,2X,A8,3X)) 

LIMUIBEG 

LIM2=LIMU(2*M-1  )*ND2 
DO  40  L=1,ND2 

WRITE  (KPTR, 45)  STORE(L),  (STORE(J) , J=LIM1 , LIM2,  ND2) 
FORMAT  (1H  ,F6.2, 10(2X, 1PE9.2)) 

LIM1=LIMU1 

LIM2=LIM2+1 

CONTINUE 

PLOT  OPTIONS 
IM1=(I-1 )»20 
DO  90  K=2,4,2 
DO  90  L=1,20 
IQ=LPRNT(L,I) 

IF  (IQ.NE.K  .AND.  IQ.NE.K-1)  GO  TO  90 
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YH=XMAX(L+IM1) 

YL=XMIN(L+IM1) 

ISTrlO 

IEND=IBEG+ND2-1 
PSYM" PS  1 

call'kplokgraph  .store,  ST0RE(  IBEG)  , 0 , 0 , 0 , 0 , 0 ) 


DO  80  J=IBEG,IEND 
C1=STORE(J) 
C2=ST0RE(J+ND2) 
ST0RE(J)=C1-C2 
ST0RE(J+ND2)=C1+C2 
80  CONTINUE 


82 

85 

90 


IST=0 
PSYM*"  PSP 

CALL'kPLOT(GRAPH, STORE, ST0RE( IBEG) , 0, 0, 0, 0, 0) 

1ST- 1 

CALL  KPLOT(GR APH , STORE , ST0RE{ IBEG+ND2 ), 0,0, 0,0,0) 
WRITE  (KPTR.82)„TITLE(I),  L 
FORMAT  (/,1H  ,A8,I3) 


IBEG=IBEG+2»ND2 

CONTINUE 


RETURN 

END 
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C 

C 


10 

C 

1600 

40 

50 

C 

100 

300 

600 
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C 

1200 

1300 

1400 

1500 

1700 


SUBROUTINE  KPLOT ( W , X , Y , NTAPE , IX , lY , NV AR , Y 1 ) 

COMMON  /PL0T1/ 

N V . NH,  NC PW , LW , VLH ( 4 ) , NXES , NDIR , 1ST , NGLV , NGLH , BS YM,  GS YM, 
PSYM,NDIM1,NDIM2,N0 

DIMENSION  W(1) ,X(1 ) ,Y(1) ,Y1(NVAR),ST0RE(70) ,0(4) ,IPX(4) ,K(3) 

EQUIVALENCE  (Q( 1 ) .XL1 ) , (0(2 ) .XH1 ) , (0(3) , YL1) , (0(4 ) , YH1 ) 
EQUIVALENCE  (ISC , K( 1 )),( JSC, K( 3 ) ) 

DATA  IPX/3,4, 1,2/ 

IF  (NH.GT.121)  NH=121 

NCPW  IS  THE  NUMBER  OF  CHARACTERS  PER  WORD 
(60  BIT  WORD  6 BIT  DISPLAY  CODE  ON  CDC) 

NCPWrIO 

LW=NH/NCPW+1 

IF  ((IST/10).GT.0)  NC0UNT=0 
NC0UNT=NC0UNT+1 
IF  (NCOUNT.EQ.10)  IST=1 
L = 1 

DO  10  1=1.4 

Q(lj=-1.0E08*(-1)«»I 

K(L)=1 

IF  (VLH(L) .EQ.VLH(L+1 ))  GO  TO  10 

K(L)=0 

Q(I)=VLH(I) 

IF  (I.EQ.2)  L=3 

IF  ( NTAPE. EQ.O)  GO  TO  1200 


SKIP  THIS  PART  IF  PLOTTING  FROM  CORE 

FLAG=0 

GO  TO  40 

IFLAG=1 

NN=0 

REWIND  NTAPE 
READ  (NTAPE)  Y1 
GO  TO  2800  ON  EOF 
IF  (EOF(NTAPE)) 2800, 100 
NN=NN-f1 

IF  (NN.LT.NDIM1)  GO  TO  50 
IF  (IFLAG.EQ.1)  GO  TO  1700 
IF  (ISC+JSC.EQ.O)  GO  TO  1710 
NN=NN+1 


IF  (ISC. EQ.O)  GO  TO  6OO 
XLUAMIN1CXL1  ,Y1(IX)) 
XH1=AMAX1(XH1,Y1(IX)) 


IF  (JSC.E(j.O)  GO  TO  200 

YL1=AMIN1(YL1,Y1(IY)1 

YH1=AMAX1(YH1,Y1UY)) 

READ  (NTAPE)  Y1 
IF  (E0F(NTAPE))1600,210 
IF  (NN-NDIM2)  300,300,1600 

RESUME  HERE 

IF  (ISC. EQ.O)  GO  TO  1400 
DO  1300  I=NDIM1,NDIM2 
XL1=AMIN1(XL1,X(I)) 
XH1=AMAX1(XH1,X(I)) 

IF  (JSC, EQ.O)  GO  TO  1700 
DO  1500  I=NDIM1,NDIM2 
YL1=AMIN1(YL1,Y(I)) 

YH1=AMAX1(YH1, Y(I)) 

IF  (ISC.EQ.1)  CALL  ADJUST(XH1 , XL1 ) 
IF  (JSC.EQ.l)  CALL  ADJUST( YH 1 , YL1 ) 
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1710 

1720 


1740 

1760 


1780 


1800 


2?00 

C 


2400 


C 

2500 

2600 

C 

2700 


C 

2800 

2900 


IF  (NDIR/10)  1720,1740,1720 

TMP=XL1 

XL1=XH1 

XH1=TMP 

IF  (NDIR-10»(NDIR/10))  1760,1780,1760 

IMP=YL1 

YL1=YH1 

YHUTMP 

J=7»(NC0UNT-1 )+1  , ^ 

IF  (J.EQ.l)  CALL  QINIT(W) 

STORE(J)=PSYM 
DO  1800  1=1,4 
IF  (NXES.EQ.O)  L=I+J  ^ 

IF  UXES.GT.O)  L=IPX(I)+J 
ST0RE(L)=Q(I) 

ST0RE(J+5)='NH-1 )/(ST0RE(J+2)-ST0RE( J+1 ) ) 
ST0RE(J+S)=(NV-1 )/(ST0RE(J+4)-ST0RE(J+3)) 
IF  (NTAPE.EQ.O)  go  to  2500 

acip  THIS  PART  IF  PLOTTING  FROM  CORE 
DO  2400  I=NDIM1,NDIM2 
IF  (NXES.EQ.O)  CALL  KPLOTC 
IF  (NXES.GT.O)  CALL  KPLOTC 
READ  (NTAPE)  Y1 
IF  (EOF(NTAPE))  2800,2700 


DO  2600  IrNDIMI . NDIM2 


(STORE(c 

I),W,Y1l 

[IX),Y1(IY)) 

(STORED 

rhw,Yii 

ilY),YlUx)) 

FROM  A 

FILE 

( STORE (c 

i),w,x(; 

(STORECt 

1),W,YC 

RESUME  HERE 

IF  { (IST-10«(IST/10)) .GT.O)  CALL  QPRINT(W, NO, NCOUNT, STORE) 
RETURN 

ERROR  MESSAGE 
WRITE  (NO, 2900) 

FORMAT(//32H  INSUFFICIENT  DATA  ON  INPUT  FILE,/, 

1H  ,28H  PLOTTING  ROUTINE  TERMINATED) 

RETURN 

END 
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SUBROUTINE  ADJUSKXHI , XL1 ) 

IF  (XH1.EQ.XL1)  XL1=0.9*XL1-10.0 
A=IFIX (100. 0+ALOG1 0 ( XH 1-XL 1 ) ) - 1 00 . 0 
XH1T=XH1*10.0»*(1 .0-A) 
XL1T=XL1*10.0«*(1 .0-A) 

IF  (XH1T  .GE.  0.0)  XH1T=IFIX(XH1T+0.9) 
XH1T=IFIX(XH1T) 

IF  (XL1T  .LE.  0.0)  XL1T=IFIX(XL1T-0.9) 
XL1T=IFIX(XL1T) 

XH1=XH1T»10.0»»(A-1 .0) 
XL1=XL1T»10.0»»(A-1 .0) 

RETURN 

END 


I 
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SUBROUTINE  QINIT( IMAGE) 

COMMON  /PL0T1/ 

1 NV , NH , NC  PW , LW , Q( 4 ) , NXES , NDIR , 1ST , NGLV , NGLH , BSYM , GSYM 
DIMENSION  IMAGE(I) 

DATA  IBLNK/10H  / 


N=LW»NV 
DO  100  1=1, N 

100  IMAGE(I)=IBLNK 
DO  101  1=1, NH 

CALL  QPLOTf IMAGE, 1,1, BSYM) 

101  CALL  QPLOT(IMAGE,I,NV,BSYM) 
DO  102  1=1, NV 

CALL  QPLOT(IMAGE,1,I,BSYM) 

102  CALL  QPLOT(IMAGE,NH,I.BSYM) 

1800  IF  (NGLV.EQ.O)  GO  TO  2000 

NGLV1=NGLV+1 

NH1=NH-1 

DO  1900  I=NGLV1,NH1,NGLV 
DO  1900  J=1,NV 

1900  CALL  QPLOT(IMAGE.I,J,GSYM) 

2000  IF  (NGLH.EQ.O)  RETURN 

NGLH1=NGLH+1 
NV1=NV-1 

DO  2100  I=NGLH1,NV1,NGLH 
DO  2100  J=1,NH 

2100  CALL  QPLOT(IMAGE, J,I,GSYM) 

RETURN 
END 
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SUBROUTINE  KPL0TC( W, IMAGE, X,Y) 

DIMENSION  W(1),  IMAGE(I) 

COMMON  /PL0T1/  NV,  NH 


IF^((J^LEloKOR!(j^GT.NH))  RETURN 
I=NV-IFIX( (Y-W(4) )*W(7)+0.5) 

IF  ((liEloKoLU.GT.NV  ) RETURN 
CALL  QPLOT  (IMAGE, J, I,W( 1)) 

RETURN 

END 


1 

[■ 
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SUBROUTINE  QPLOT(IMAGE, J, I.SYM) 

COMMON  /PL0T1/  NV,  NH,  NCPW,  LW 
DIMENSION  IMAGE(I) 


101 

102 

103 


II=J/NCPW 

L=J-NCPW»II 


11=11+1 

IF  (L)  101,101,102 
L=NCPW 
11=11-1 


IW“II+(I-1 )*LW 

aLL  PLACE(IMAGE(IW),L,SYM, 

RETURN 

END 


1) 
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SUBROUTINE  PUCE(  A,  N,  B.M) 

C THE  MTH  CHAR  OF  B REPLACES 

C THE  NTH  CHARACTER  OF  A 

C CHAR  POSITIONS  ARE  1 TO  10  FROM  LEFT  TO  RIGHT 

COMMON/ I NOU/KIN, KOUT 
INTEGER  A,  B,  BX,  BY 
DATA  MASK/77B/ 

C CHECK  FOR  VALID  ARGUMENTS 

IF  (N.GT.10  .OR.  M.GT.10)  GO  TO  900 

IF  (N.LT.1  .OR.  M.LT.1)  GO  TO  900 

C NULL  ALL  BUT  THE  MTH  CHAR  OF  B,  PUT  IT  IN  BX 

C NULL  THE  NTH  CHAR  OF  A 

NSHFT=60-6»N 
MSHFT=60-6»M 

MASKBY  = SHIFT(MASK,NSHFT) 

MASKB  = SHIFT ( MASK, MSHFT) 

MASKA  = COMPL(MASKBY) 

A = AND(A.MASKA) 

BX  = AND(B, MASKB) 

C SHIFT  THE  MTH  CHAR  OF  BX  TO  THE  NTH  POSITION 

C AND  PUT  IT  IN  BY 

MNSHFT=6»(M-N) 

BY  = JHIFT(BX,MNSHFT) 

BY  = AND( BY, MASKBY) 

C COMBINE  A AND  BY 

A = OR(A,BY) 

RETURN 

C N OR  M OUT  OF  BOUNDS 

900  WRITE  (KOUT, 1900) 

1900  FORMAT  (1H  .20HERROR  IN  SUBR.  PLACE,/, 

1 24H  N OR  M IS  OUT  OF  BOUNDS) 

CALL  EXIT 
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110 


120 


130 

150 

102 

103 

105 


SUBROUTINE  QPRINT( IMAGE, NO, NCOUNT, STORE) 

DIMENSION  IMAGE(I),  STORE(I) 

COMMON  /PL0T1/  NV,  NH,  NCPW,  LW 

CALL  PAGEFD(N0, 1) 

DO  110  1=1, NCOUNT 

WRI^E^InoI^W)  STORECIB+1  ) ,ST0RE(IB)  ,ST0RE(IB)  ,ST0RE(IB+2) 

NCANT=NV-NCOUNT 

IA=1 

DO  150  1=1, NV 
IB=I»LW 

IF  (I. GT. NCOUNT)  GO  TO  120 

IBASE=(I-1)«7+1  , , 

WRITE  (NO, 103)  ST0RE(IBASE),ST0RE(IBASE+4),(IMAGE(J),J=IA,IB) 
GO  TO  150 

IF  (I.GT.NCANT)  GO  TO  I30 
WRITE  (NO, 105)  (IMAGE(J),J=IA,IB) 

GO  TO  150 

IBASE=(I-1-NCANT)*7+1  . , . 

WRITE  (NO, 103)  ST0RE(IBASE),ST0RE(IBASE+3),(IMAGE(J),J=IA,IB) 
IA=IA+LW 

FORMATTiH  , 11X,1PE10.3, 1X.A1.57X.A1, 1X,1PE10.3) 

F0RMAT(1H  ,A1. 1PE9.2, IX, 12A10,A1) 

FORMATUH  , 1 1X,  12A10, A1 ) 

RETURN 

END 
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LIBRARY  FOR  OPTIMAL  CONTROL  MODEL 


SUBROUTINE  KFLTR ( N , M . A . yT . Q . R , X » G . Z ) 

SOLVES  X=AfX-XH  (HXH  +RThX)A  +Q  , _ , , , 

DIMENSION  A(1),  HT(1),  Q(l),  R(1),  X(1),  G(1),  Z( 1 ) , TR(80) 
COMMON  /MAIN1/  NDIM,  NDIM1,  F( 1 ) 

COMMON  /INOU/  KIN,  KOUT 

NN=N»NDIM 

NT=1 

11=0 

DO  1 I = 1,N 

IF  (2»»NT.GT.N)  GO  TO  2 

NT=NT+1 

DO  4 J=1,M 

DO  3 1=1, N 

K=II+I 

G(K)=HT(K)/R(J) 

II=II-*-NDIM  , 

CALL  EQUATE(X,A,N,N) 

CALL  MAT2(N,M,G,HT,Z) 

DO  15  L=1,NT 

CALL  MAT3A(N,N,X,Z,F) 

DO  10  1=1, N 

DO  10  J=I,NN,NDIM 

Z(J)=Z(J)+F(i) 

cal^'mmu^( F, F, N, N, N, X) 

CALL  F ACTOR (N,Z,F, MR) 

IF  (MR.EQ.-1)  CALL  EXIT 
CALL  GMINV(N,N,F,Z,MR,0) 

WRITE  (KOUT, 16)  N,  N,  MR 

FORMAT  (25HOOBSERVABILITY  MATRIX  IS  , 12, 1HX, 12, 9H  OF  RANK  ,12) 
CALL  MMUL(X,Z,N,N,N,F) 

CALL  MAT2(N,N,F,F,X) 

ENTRY  KFLTR 1 
TOL=1 .E-5 
ADV=T0L»1 .E-7 
NN=N»NDIM 
L=0 

TOLUTOL/10. 

MAnT=30 
DO  19  1=1. N 
TR(I)=-1 .0 
CONTINUE 

DO  40  IT=1.MAXIT 
CALL  MMUL( X , HT, N, N, M, F) 

CALL  MAT2A(N,M,HT,F,G) 

11=1 

DO  20  1=1, M 
G(II)=G(II)+R(I) 

II=II+NDIM1 

CALL  GMINV(M,M,G.Z,MR,0) 

CALL  MMUL(F,i,N,M,M,G) 

CALL  MAT5(G,HT,N,M,N,Z) 

11=1 

DO  18  1=1, N 

DO  17  J=I,NN,NDIM 

Z(J)=-Z(J)^ 

Z(II)=Z(II)+1 .0 
II=II+NDIM1 
IF  (L.EQ.N)  RETURN 
CALL  MMUL(A,Z,N,N,N,X) 

CALL  MMUL(A,G,N,N,M,F) 

11=0 


1 
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DO  22  J=1,M 
DO  21  1=1, N 
K=II+I 

21  G(K)=F(K)*R(J) 

22  II=II+NDIM 

CALL  MAT2(N,M,G,F,Z) 

CALL  MADD1(N,N,Z,Q,Z,1.0) 

CALL  TRANS1(N,X,X) 

CALL  DLINEQ(N,X,Z,X,T0L1) 

L=0 

C1=0.0 

11=1 

DO  25  1=1, N 

IF  (ABS(X(II)-TR(I)).LT.(ADV+TOL»X(II)))  L=L+1 
TR(I)=X(II) 

II=II+NDIM1 
25  C1=C1+TR(I) 

IF  (ABS(Cl)  .GT.  1.0E20)  GO  TO  50 
IF  (L.NE.N)  GO  TO  40 
WRITE  (KOUT.27)  IT 

27  format} 17H0RICCATI  SOLN  IN  ,I2,11H  ITERATIONS) 

CALL  GMINV(N,N,Z.F,MR,0) 

IF  (MR.NE.N)  WRITE  (KOUT.35)  MR 

35  F0RMAT(26H0RICCATI  SOLN  IS  PSD— RANKI3) 

40  CONTINUE 

WRITE  (K0UT.45)  MAXIT 

45  FORMAT} 26H0RICCATI  NON-CONVERGENT  INI2.11H  ITERATIONS) 

CALL  EXIT 

50  WRITE  (KOUT.55)  IT,  NT 

55  F0RMAT}29H0RICCATI  BLOW  UP  AT  ITERATIONI2, 12H  INITIAL  N=I3) 

CALL  EHT 
RETURN 
END 
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SUBROUTINE  MLINEQ(N, A, C.X.TOLI ) 
ANSWER  RETURNED  IN  C AND  X 
DIMENSION  A(1),C(1),X(1) 

COMMON  /MAIN1/  NDIM,  NDIM1,  F( 1 ) 
COMMON  /INOU/  KIN,  KOUT 

DT=.5 

DT1=0. 

NN=N*NDIM 

DO  5 II=1,NN,NDIM1 

DT1=DT1+ABS(A(II)) 

DTUDT1/N 

IF  (DTI. GT. 4.0)  DT=DT»4.0/DT1 
11=1 

DO  20  1=1, N 

DO  15  JJ=1,NN,NDIM 

X(JJ)=DT»A(JJ) 

X(II)=X(II)-.5 

II=II+NDIM1 

CALL  GMINV(N,N.X,F,MR, 1) 

IF  (MR.EQ.N)  GO  TO  21 
IT=0 

DO  18  I=1,NN,NDIM1 
C(I)=1.E25 
GO  TO  95 

CALL  MMUL(C,F,N,N,N,X) 
INITIALIZATION  OF  X,F 
1=1 

DO  40  II=1,NN,NDIM 
J=II 

IF  (I.EQ.1)  GO  TO  30 

J=J+1 

ID=J 


DO  35  JJ=II,NN,NDIM 
C(J)=DT»DOT(N,F(II) ,X(JJ)) 

J=J+1 

F(ID)=F(ID)+1 .0 
1=1+1 
GO  TO  50 
ENTRY  DLINEQ 
NN=N«NDIM 

CALL  EQUATE(F,A,N,N) 

TOL=TOL1 

ITT=31 

IF  (TOL1  .GT.  0.99)  ITT=IFIX(TOL) 

IF  (TOL1  .GT.  0.99)  TOL=0.0 

ADV=TOL»1 .E-7 

DO  90  IT=1,ITT 

NEZ=0 

SIZE=0.0 

CALL  MMUL(C,F,N,N,N,X) 

ii=i 

J=1 

GO  TO  70 
J=II 

DO  65  JJ=I,II,NDIM 
C(J)=C(JJ) 

J=J+1 

ID=J 

DT1=C(J) 


DO  75  JJ=II.NN. 
C(J)  = C(J)+D0T(N 
J=J+1 


NDIM 

,F(II),X(JJ)) 
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80 


90 

95 

100 

150 


,LT. ( ADV+T0L«ABS( C( ID) ) ) ) NEZ=NEZ+1 


J=J-1 

DO  80  JJ=II,J 
X(JJ)=F(JJ)  , 

IF  (ABS(C(ID)-DTI) 

1=1+1 
II=II+NDIM 
SIZE=SIZE+DT1 
IF  (I.LE.N)  GO  TO  60^^^ 

IF  1aBS(IiZe!  °Gt!  1 Jei8)  GO  TO  95 

CALL  MMUL(X,X.N,N,N,F) 

IF  (IT.EQ.ITT)  GO  TO  150 
CONTINUE  , 

FoSl‘(33H6uS'Egi  algorithm  HOH-COHVERGE»I,I3,10HITEBATIOHS) 
CALL  EQUATE(X,C,N,N) 

RETURN 

END 
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SUBROUTINE  DINTEG(N . A, C*S, NT)  . 

C S=SUM  1=0  TO  NT-1  OF  A I*C»AT  I 

DIMENSION  A(1),  C( 1 ) , S(1) 

COMMON  , , 

1 /MAIN1/  NDIM,  NDIM1,  COMI(I) 

2 /MAIN2/  C0M2(1) 

KsNT 

1MX=(N-1 )«NDIM+N 

DO  10  1=1, N 

DO  10  JrI,NMX,NDIM 

^J)=0.0 

C0M2(J)=C(J) 

10  CONTINUE 

IFLAG=0 

100  IP=M0D(K,2) 

K=K/2 

IF  (IP.EQ.O)  GO  TO  110 

IF  (IFLAG.EQ.O)  CALL  EQUATE{COM1, A. N.N) 

IF  UFLAG.EQ.I)  CALL  MMUL( A, C,N, N, N, COM1 ) 
IFUG=1 

DO  14  1=1, N 

DO  14  J=I,NMX,NDIM 

14  CONTINUE 

IF  (K.EQ.O)  GO  TO  200 

CALL  MMUL(A,C0M2,N,N,N.C0M1) 

CALL  MAT2(N,N.COM1,A,C0M2)  ^ 

110  CALL  MMUL(A,c6M2,N,N,N,COM1) 

CALL  MAT6S(N,N,C0M1,A.C0M2) 

CALL  MMUL(A,A,N,N,N,C0M1) 

CALL  EQUATEU,COM1,N,N) 

GO  TO  100 


200 


RETURN 

END 
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SUBROUTINE  DSCRT( N , A , DEL , EA , EAINT , NT ) 

DIMENSION  A( 1 ) ,EA( 1 ) ,EAINT( 1 ) ,C0EF(30) 

C SETS  EA=EXP(A«DEL) .EAINTrlNTEGRAL  EA  0 TO  DEL 

C OMMON/MAIN 1 /NDIM , NDIM 1 
NN=N*NDIM 
NTM1=NT-1 
COEF(NT)=1 . 

DO  10  I=1,NTM1 
II=NT-I 

10  COEF(II)=DEL«COEF(II+1)/FLOAT(I) 

C NT  MUST  BE  AT  LEAST  3 

11=1 

DO  30  1=1 ,N 
DO  20  J=I,NN,NDIM 
20  EAINT(J)=A(J)*COEF(1 ) 

EAINT(II)=EAINT(II)+C0EF(2) 

30  II=I1+NDIM1 

DO  60  L=3,NT 
T1sCOBF(L) 

CALL  MMUL(A.EAINT,N,N,N,EA) 

ifTl.eq.nt)go  to  70 
11=1 

DO  60  1=1, N 
DO  50  J=I,NN,NDIM 
50  EAINT(J)=EA(J) 

EAINT(II)=EAINT(II)+T1 
60  II=II+NDIM1 

70  DO  80  II=1,NN,NDIM1 

EA(II)=EA(II)+T1 
80  CONTINUE 

RETURN 


END 
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1000 

IE 

42 


43 

50 

C 

55 


60 


SUBROUTINE  FACTOR (N, A, S, MR) 

A=ST*S 

DIMENSION  A(1),  S(1),  IP(80) 

COMMON  /MAIN1/  NDIM,  NDIM1 
COMMON  /INOU/  KIN,  KOUT 
MR=0 

NN=N*NDIM 
T0L=1 .OE-06 

CALL  MSCALE(S,A,N,N,1.0) 

1=1 

DO  50  ID=1,NN,NDIM1 

NEND=ID-I 

R=Q.O 

IP(I)=I 

K=I 

DO  20  KK=ID,NN,NDIM1 

t=sTkk)-dot2(nend,s(k),s(P) 

IF  (ABS(T)  .LE.  ('{'0L»(AbS(S(KK))+T0L)))  GO  TO  20 
IF  IkESiT)  .LE.  ABS(R))  GO  TO  20 
IP(I)  = K 
R=T 

K=K+1  , ^ ^ 

IF  (I  .EQ.  IP(D)  GO  TO  25 
K“IP(I) 

CALL  SWAP(S(I),S(K),N,NDIM) 

K=ID+NDIM«(K-I) 

NP=N-I+1 

CALL  SWAP(S(ID),S(K),NP,1) 

NP=ID 

IF  (R)  30,  35,  40 
MR=-1 

WRITE  (KOUT, 1000) 

FORMAT  (1H  ,33HTRIED  TO  FACTOR  INDEFINITE  MATRIX) 

RETURN 

NP=NEND+N 

K=NEND+1 

DO  42  KK=K,NP 

S(KK)=0.0 

IF  (R  .EQ.  0.0)  GO  TO  50 
T=SQRT(R) 

S(ID)  = T 
MR=MR+1 

IF  (I  .EQ.  N)  GO  TO  55 

ID1=ID+1 

K=I+1 

NP=NEND+N 
DO  43  KK=ID1,NP 

S(KK)=(S(KK)-D0T2(NEND,S(I),S(K)))/T 

K=K+1 

1=1+1 


UNSWAP  ROWS 
DO  60  1=1, N 
K=N-I+1 
KK=IP(K) 

IF”(KK  .EQ.  K)  GO  TO  60 

CALL  SWAP(S(K) ,S(KK) ,MR,NDIM) 

CONTINUE 

CALL  TRANS1(N,S,S) 

RETURN 

END 


1 

I 
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SUBROUTINE  GMINV ( NR , NC , A , U , MR , MT ) 

DIMENSION  A(1),U(1),S(30) 

C0MM0N/MAIN1/  NDIM.NDIMI 

COMMON/INOU/  KIN,  KOUT 

T0L=1 .E-12 

MR=NC 

NRMUNR-1 

T0L1=1 .E-20 

JJ=1 

DO  100  J=1,NC 
FAC=DOT(NR,A(JJ) ,A(JJ)) 

JMlrJ-1 
JRM=JJ+NRM1 
JCM=JJ+JM1 
DO  20  I=JJ,JCM 
U(I)=0. 

U(JCM)=1.0 
IF(J.EQ.I)  GO  TO  54 
KK=1 

Da"39  Ksi.JMi 
IF(S(K).EQ.1.0)  GO  TO  30 
TEMP=-DOT(NR,A( JJ) ,A(KK) ) 

CALL  VADD(K,TEMP,U(JJ),U(KK)) 
KK=KK+NDIM 
DO  50  L=1,2 
KK=  1 

DO'50  K=1,JM1 

IF(S(K) .EQ.O.)  GO  TO  50 

TEMP=-DOTCNR,A(JJ>,A(KK)) 

CALL  VADD{NR.TEMP,A(JJ),A(KK)) 

CALL  VADD(K,TEMP,U(JJ) ,U(KK)) 

KK=KK+NDIM 

TOL1=TOL«FAC 

FAC=DOT(NR,A(JJ),A(JJ)) 

IF(FAC.GT.TOLI)  GO  TO  70 

DO  55  I=JJ,JRM 

A(I)=0. 

S(J)=0. 

KK“  1 

do" 65  K=1,JM1 

IF(S(K) .e6.0.)  go  to  65 

TEMP=-D0T?K, U(KK) ,U( JJ)) 

CALL  VADD(NR,TEMP,A(JJ),A(KK)) 
KK=KK+NDIM 

FAC=DOT(J,U(JJ) ,U(JJ)) 

MR=MR-1 
GO  TO  75 
S(J)=1.0 
KK"1 

DO'72  Ksl.JMI 
IF(S(K).EQ.1.)  GO  TO  72 
TEMP=-DOT(NR,A(JJ),A(KK)) 

CALL  VADD(K,TEMP,U{JJ) ,U(KK)) 

KK=KK+NDIM 

FACs T./SQRT( FAC) 

DO  80  IsJJ.JRM 
A(I)sA(I)*FAC 
DO  85  IsJJ, JCM 
U(I)sU(I)»FAC 
JJsJJ+NDIM 

IF(MR.EQ.NR.OR.MR.EQ.NC)  GO  TO  120 
IF(MT.NE.O)  WRITE  (KOUT, 110)  NR,NC,MR 
F0RMAT(I3,1HX,I2,8H  M RANK, 12) 

NENDsNC*NDIM 

JJsl 


,TEMP,U(JJ) ,U(KK)) 
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DO  135  J = 1,NC 
DO  125  1=1, NR 
II=I-J 
S(I)=0. 

DO  125  KK=JJ,NEND,NDIM 
125  S(I)=S(I)+A(II+KK)«U(KK) 

II=J 

DO  no  1=1,  NR 
U(II)=S(I) 

130  II=II+NDIM 

135  JJ=JJ+NDIM1 

RETURN 
END 
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SUBROUTINE  MAT2(N1 .N2,X, Y,Z) 
Z=XY  X,Y=N1»N2,Z=Z 
Z AND  Y CAN  BE  EQUIVALENT 
DIMENSION  X(1).Y(1),Z(1) 
C0MM0N/MAIN1 /NDIM , NDIM1 
NN2=N2*NDIM 
11=1 

DO  10  1=1, N1 

ij=n 

DO  5 J=I,N1 

Z(IJ)=D0T2(NN2,X(I),Y(J)) 

IJ=IJ+NDIM 

J=II 

LJ=J 

IJ=IJ-NDIM 
IF(IJ.LT.I)  GO  TO  10 
J=J-1 

Z(IJ)=Z(J) 

GO  TO  3 
II=II>NDIM1 
RETURN 
END 


si 

j 


oo 
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SUBROUTINE  MAT2A(N1 ,N2,X,Y,Z) 
Z=XT»Y  X,Y  ARE  N1XN2,  Z=ZT 
Z AND  Y CAN  BE  EQUIVALENT 
DIMENSION  X(1).Y(1),Z(1) 
COMMON/MAIN 1/NDIM 
NN2=N2*NDIM 
1=0 

DO  10  II=1,NN2,NDIM 
J=II+I 


U=J-1 

DO  5 JJ=II,NN2.NDIM 
Z(J)=D0T(N1,X(1I),Y(JJ)) 
J=J+1 
1=1+1 
JJ=I 

DO  10  J=II,IJ 
^J)  = Z(JJ) 

JJ=JJ+NDIM 

RETURN 

END 
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SUBROUTINE  MAT3(N1 , N2, X, Y , 
Z=XYXT  Y=YT  IS  N2XN2^ 
DIMENSION  X(1),Y(1),Z(1) 
C0MM0N/MAIN1/  NDIM 
CALL  MMUL(X,Y,N1,N2,N2,Z) 
NN2=N2»NDIM 
DO  5 r=1,N1 
IM1=I-1 
II=IM1»NDIM 
JJ=I+II 

DO  3 J=I,N1  ^ 

Z(JJ)  = D0T2(NN2,Z(J)  ,X(D) 

JJ=JJ+NDIM 

JJ=I 

K=II-*-1 

KK=II+IM1 

DO  5 J=K,KK 

Z(JJ)=Z(J) 

JJ=JJ+NDIM 

CONTINUE 

RETURN 

END 


Z) 

X IS  N1XN2 
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SUBROUTINE  MAT5(X,Y,N1 ,N2,N3,Z) 
DIMENSION  X(1) ,Y(1),Z(1) 
C0MM0N/MAIN1/NDIM , NDIM1 

Z=X«YT  X IS  N1«N2,  Y IS  N3*N2 
JJ=1 

DO  3 I=1,N3 
C1=Y(I) 

CALL  VSCALE(Z(JJ),X,N1,C1) 

JJ=JJ+NDIM 

KrNDIMI 

GO  TO  4 

ENTRY  MAT5S 

K = 1 

NN2=N2*NDIM+1 


NN3=N3*NDIM 

IF(K.EQ.NN2 

I=K 


EQ.NN2)  RETURN 


DO  6 JJ=1,NN3,NDIM 
C1“Y(I) 

IF(C1.NE.O.)  CALL  VADD(N1 , Cl , Z( JJ) ,X(K) ) 
1=1+1 
K=K+NDIM 
GO  TO  7 
END 
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SUBROUTINE  MAT5A(X, Y, N1 . N2, N3, Z) 

Z=XT«y  X=N2*N1,  Y=N2»N3 
DIMENSION  X(1),Y(1),Z(1) 

COMMON/MAIN 1/NDIM 

N1MUN1-1 

NN3=N3*NDIM 

DO  1 I=1,NN3,NDIM 

II=I+N1M1 

DO  1 J=I,II 

Z(J)=0.0 

ENTRY  MAT5AS 

m3=N3*NDIM 

DO  10  K=1,N2 

KK=K 

DO  8 1=1, N1 

IF(CKNi.O.O)  CALL  VADD1  (NN3,  Cl  ,Z(  I) , Y(K) ) 

KK=KK+NDIM 

CONTINUE 

RETURN 

END 


I 


Report  No . 3M6^» 


Bolt  Beranek  and  Newman  Inc. 


C 

1 

C 

6 

8 

10 


SUBROUTINE  MAT6(N1 , N2, X, Y, Z) 

Z=X»Y  . WHERE  X=N1»N2,  Y=N1»N2,  Z=Z  =N1»N1 
DIMENSION  X(i;,  Y(1),  Z( 1 ) 

COMMON  /MAIN1/  NDIM,  NDIM1 

NN1=N1»NDIM 

DO  1 1=1, N1 

DO  1 J=I,NN1,NDIM 

Z(J)=0.0 

CONTINUE 

ENTRY  MAT6S 

Z=Z+X*Y 

NN2=N2*NDIM 

NN1=N1«NDIM 

DO  6 K=1,NN2,NDIM 

KK=K-1 

J=1 

DO  6 1=1, N1 
C1=Y(I+KK) 

IF  (C1.NE.0.0)  CALL  VADD(I ,C1 , Z( J) ,X(K) ) 

J=J+NDIM 

CONTINUE 

IF  (N1.EQ.1)  RETURN 

NN2=NDIM1+1 

DO  10  K=NN2,NN1,NDIM1 

I=K 

J=K 

I=I-t 

J=J-NDIM^ 

z(j)=zfi) 

IF  (J.GT.NDIM)  GO  TO  8 

CONTINUE 

RETURN 

END 


r 

[ 

t 

I 

I 

I 

i 


[ 

f 

t 

[ 

t 
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SUBROUTINE  MMUL(X,Y,N1 ,N2,N3, 

DIMENSION  X(1),Y(1) ,Z(1) 

C0MM0N/MAIN1/NDIM 

N1M1=N1-1 

NN3=N3*NDIM 

DO  1 I=1,NN3,NDIM 

II=I+N1M1 

DO  1 J=I,II 

Z(J)=0.0 

ENTRY  MMULS 

NN^=N^*NDIM 


Z) 


KK=0 

DO  10  K=1.N2 
DO  8 1=1. N1 
I+KiC) 

IF(CI.NE.O.O)  CALL  VADD1 (NN3.C1 ,Z(I) ,Y(K) ) 

CONTINUE 

KK=KK+NDIM 

RETURN 

□ID 
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SUBROUTINE  MADD1 (NR, NC, X, Y, Z,C1) 
C Z=X+C1»Y 

DIMENSION  X(1),Y(1),Z(1) 
C0MM0N/MAIN1/NDIM 
NN=NC»NDIM 
IF(C1.EQ.1.0)G0  TO  1 
IF(C1.EQ.-1  JGO  TO  2 
DO  5 1=1, NR 
DO  5 J=I,NN,NDIM 
5 Z(J)=X(J)+C1»Y(J) 

RETURN 

1 DO  10  1=1, NR 

DO  10  J=I,NN,NDIM 
10  Z(J)  = X(J)+Y(J) 

RETURN 

2 DO  15  1=1, NR 

DO  15  J=I,NN,NDIM 
15  Z(J)=X(J)-Y(J) 

RETURN 

END 


oo 
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SUBROUTINE  DIAG2(N,A,B,C1 ,C2) 

A = C1«B  + C2»I 

A,B  ARE  N«N  MATRICES;  I IS  N«N  IDENTITY  MATRIX 
DIMENSION  A(1),  B(1) 

COMMON  /MAIN1/  NDIM,  NDIM1 
NN=N»NDIM 
NM1=N-1 
11=1 


IF 


(Cl  .EQ.  1 
DO  5 J=1,NN,N 
K=J+NM1 

if  X 

A(I)=C1»B(p 
A(II)=A(II)+C2 
II=II+NDIM1 
RETURN 

DO  7 J=1,NN,NDIM 
K=J>NM1 
DO  6 I=J,K 


0)  GO  TO 
NDIM 


10 


V/  f r 


D+C2 
Ii=Ii+NDIM1 
RETURN 
END 


’ f 


I 
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SUBROUTINE  IDENT(N,A,C1 ) 

DIMENSION  A(1) 

C0MM0N/MAIN1/NDIM,NDIM1 

NN=N«NDIM 

11=1 

DO  1 1=1, N 

DO  2 J=I,NN,NDIM 

A(J)=0.0 

A(II)=C1 

II=II+NDIM1 

RETURN 

EXD 
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SUBROUTINE  EQUATE( A, B, NR, NC) 
A?B 

MATRIX  EQUATE 
DIMENSION  A(1).  B(1) 

CALL  MSCALE(A,B,NR,NC, 1.0) 

RETURN 

END 


I 
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SUBROUTINE  MSCALE(A, B, NR, NC, Cl ) 
A=C1»B 

A AND  B MAY  BE  EQUIVALENT 
DIMENSION  A(1),  B( 1 ) 

COMMON  /MAIN1/  NDIM 


NN: 

:NC»NDIM 

1.0) 

IF 

(Cl 

.EQ. 

GO 

TO 

10 

IF 

(Cl 

.EQ. 

0.0) 

GO 

TO 

20 

IF 

lc^ 

.EQ. 

-1.) 

GO 

TO 

30 

DO 

5 1= 

I, NR 

DO 

5 J = 

:I,NN, 

,NDIM 

RETURN 

DO  15  1=1, NR 
DO  15  J=I,NN,NDIM 

a(j)=b(j) 

RETURN 

DO  25  1=1, NR 

DO  25  J=I,NN,NDIM 

A(J)=0.0 

RETURN 

DO  35  1=1, NR 
DO  35  J=I,NN,NDIM 
A(J)=-B(J) 

RETURN 

END 
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SUBROUTINE  TRANS1(N,A,AT) 

SETS  AT  = ATRANSPOSE. 

A AND  AT  CAN  BE  EQUIVALENT;  BOTH  ARE  SQUARE. 
DIMENSION  A(1),AT(1) 

COMMON/MAIN 1 /NDIM , NDIM 1 

NNsN*NDIM 

DO  1 I=1,NN,NDIM1 

U=I 

DO  1 JI=I,NN,NDIM 
TEMP=A(IJ) 

AT(IJ)=A(JI) 

AT(JI)=TEMP 

IJ=IJ+1 

CONTINUE 

RETURN 

END 


ooo 
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SUBROUTINE  TRANSP(M1 ,N2,X,XP0SE) 
COMPUTE  MATRIX  TRANSPOSE 

ypQcp.v  ' 

X ANd'xPOSE  CpNOT  BE  EQUIVALENT 
DIMENSION  X(1)  .XPOSEd) 

COMMON  /MAINl/NDIM 

ENTRY  TRANS2 

K=1-NDIM 

DO  10  1=1, N1 

K=K+NDIM 

U=I 

DO  10  J=1,N2 
XP0SE(JI)=X(IJ) 

IJ=IJ+NDIM 

JI=JI+1 

CONTINUE 

RETURN 

END 


DOUBLE  PRECISION  DDT1,  DBLE 
DIMENSION  A(1),B(1) 
DDTlsO.ODO 

IF  (NR  .LE.  0)  GO  TO  2 
DO  1 1=1. NR 

1 DDT1=DDTUDBLE(A(I)»B(D) 

2 DOTrDDTI 
RETURN 
END 


i 
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SF'fiS-TE.  0)  00  10  E 

dStLdDT|K(ACI)*B(I)) 

D0T2=DDT2 

RETURN 

END 
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FUNCTION  D0T3(N,A,B) 

DOUBLE  PRECISION  DDT3,  DBLE 
DIMENSION  A(1),B(1) 

COMMON  /MAIN1/  NDIM 
DDT3:0.0D0 

IF  (N  ,LE.  0)  GO  TO  2 
11=1 

DO  1 I=1.N 

DDT3=DDT^DBLE(A(II)*B(I)) 

Ilrfl-t-NDIM 

DOT3=DDT3 

RETURN 

QiD 
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SUBROUTINE  VADD(N,C1 ,A,B) 
DIMENSION  A(1) ,B(1) 

DO  1 1=1, N 
A(I)=A(I)+C1*B(I) 

RETURN 

END 
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SUBROUTINE  VADD1(NN,C1,A,B) 
DIMENSION  A(1),B(1) 
COMMON/MAIN 1/NDIM 

RETURN 
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SUBROUTINE  VSCALE(X, Y,N,C1) 
DIMENSION  X(1),Y(1) 

L=0 

IF(Cl.EQ.l.O)  GO  TO  5 
IF(CI.EQ.O.O)  GO  TO  8 
IF(C1.EQ.-1 .)  GO  TO  13 
L=L+1 

X(L)  = Cin(L) 

IF(L.LT.N)  GO  TO  1 

RETURN 

L=L+1 

X(L)=Y(L) 

IFCL.lt. N)  GO  TO  5 

RETURN 

L=L+1 

IFl‘L'iLT?N)  GO  TO  8 

RETURN 

L=L+1 

SS'I^Lt!!)!  go  TO  13 

RETURN 

END 
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SUBROUTINE  SWAP(A,B,N, INC) 


SNAPS  A ROW,  COLUMN, 
INC=1  ROW 

INC=NDIM  COLUMN 

INC=NDIM+1  DIAGONAL 
DIMENSION  A(1),  B(  1 ) 
NNsN*INC 


6r  diagonal  of  TWO  MATRICES 


1=^ 
IF  (I 


TQ1P=A(I) 

A(I)=BU) 

B(I)=TEMP 

I=I+INC 

GO  TO  1 

END 


.GT.  NN)  RETURN 


100  - 
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« 

I 


SUBROUTINE  VMAT1(A,X,N1 ,N2, Y) 

z y =Ax 

DIMENSION  A(1),X(1),Y(1) 

COMMON/MAIN 1/NDIM 

DO  1 1=1, N1 

Y(I)=0.0 

II=I 

DO  1 J=1,N2 
Y(I)=Y(i5+A(II)«X(J) 

1 II=II+NDIM 

RETURN 
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SUBROUTINE  VMAT2(Z,A,X,N1 ,N2,Y) 
Y=Z+AX 

DIMENSION  A(1),X(1),Z(1),Y(1) 
C0MM0N/MAIN1/NDIM 
DO  1 1=1, N1 
Y(D  = Z(l) 

y?i)=y(i)Ja(ii)»x(j) 

II=II+NDIM 

RETURN 

DID 


- 102  - 
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FUNCTION  XGAIN(TH,XM,XS) 

DIMENSION  A(5) 

DATA  A/. 225B368,-. 2521287, 1 .259695,-1 .287822, .9^06461/ 

IF  (TH.GT.O.)  GO  TO  2 

XDAINrI.O 

RETURN 

2 Y=XM 

NS=2 

IF(XS.LT.1.0E-10)XS=1 .0E-10 

IFU.EQ.O.)  NS=1 

ANS=0. 

RMS=XS**2+XM*«2 
DO  1 I=1,NS 
Z=.707*(TH+Y)/XS 
TEMP=EXP(-Z»»2) 

P=X» ( i J ( a{ If 2xJa( 4 ) ?^x1a( 3 ) ) »X+A( 2 ) ) «X+A ( 1 ) ) • 1 . 1 28379 
ERFsl .-P»TEMP 
IF  (Z.LT.O.)  ERF=-ERF 

ANS=ANS+(RMS+TH*Y)»(1.-ERF)-XS»Y»TEMP».7975 
1 Y=-Y  , ^ 

XGAIN=ANS/RMS/FLOAT(NS) 

IF(XGAIN.LT.1.E-6)  XGAIN=1.E-6 

RETURN 

END 
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100 

1000 


SUBROUTINE  OUTPUT (M,N,X. LA)  , 

MATRIX  OUTPUT  ROUTINE  (WITH  HEADER) 

DIMENSION  X(1) 

COMMON  /MAIN1/  NDIM 

COMMON  /INOU/  KIN,  KOUT,  KPTR 

DATA  KB/1H  / 

ENTRY  MATTYP 
IPIUKOUT 
GO  TO  100 
ENTRY  MATPRT 
CONTINUE 
ENTRY  OUTMAT 
KFILsKPTR 

IF  (U  .NE.  KB)  WRITE  (KFIL.1000)  LA 
FORMAT(/,1H  ,A10,8H  MATRIX  ) 

WRITE  THE  MATRIX 
10=3 

CALL  MATIO(X,M,N,IO) 

RETURN 

END 
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100 
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L 


SUBROUTINE  VECOUKX, N, LET) 

VECTOR  OUTPUT  ROUTINE  (WITH  HEADER) 


DIMENSION  X(1) 

COMMON  /MAIN1/  NDIM 

COMMON  /INOU/  KIN,  KOUT,  KPTR 

DATA  KB  /1H  / 


ENTRY  VECTYP 
KFILsKOUT 
GO  TO  100 
ENTRY  VECPRT 

KFIL=KPTR  , ^ 

IF  (LET  .NE.  KB)  WRITE  (KFIL.1000)  LET 
FORMAT  (/,1H  ,A10,8H  VECTOR  ) 


WRITE  THE  VECTOR 
10=3 

CALL  VECTIO(X,N,IO) 


RETURN 

END 
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SUBROUTINE  MATI0(X, NR, NC, 10) 

BATCH  ORIENTED  MATRIX  I/O 
10=1  INPUT  ONLY 

10=2  INPUT  AND  OUTPUT 

10=3  OUTPUT  ONLY 
10=4  PUNCH 

DIMENSION  X(1) 

COMMON  /MAIN1/  NDIM 

COMMON  /INOU/  KIN,  KOUT,  KPTR,  KPUNCH 

JEND=NC»NDIM 

GO  TO  (5,5,20,40)  10 

ceeeeeeaiHPur 

5 DO  10  1=1, NR 

READ  (KInIiOOO)  {X(IJ),  IJ=I, JEND,NDIM) 

10  CONTINUE 

IF  (10  .EQ.  1)  RETURN 

CeeeeeeaoUTPUT 

20  DO  30  1=1, NR 

WRITE  (KO(lT,2000)  (X(IJ),  IJ=I, JEND.NDIM) 
30  CONTINUE 

RETURN 

caaaeeaepuNCH 

WITE  ?kMcH,3000)  (X(IJ)  ,IJ=I,  JEND.NDIM) 
50  CONTINUE 

RETURN 

1000  FORMAT  (8E10.0) 

2000  FORMAT  (1H  .1P10E13.3) 

3000  FORMAT  (1P8E10.3) 

END 


- 106 
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SUBROUTINE  VECTIO(X, N, 10) 

BATCH  ORIENTED  VECTOR  I/O 
10=1  INPUT  ONLY 

10=2  INPUT  AND  OUTPUT 

10=3  OUTPUT  ONLY 
10=4  PUNCH 

DIMENSION  X(1) 

COMMON  /INOU/  KIN,  KOUT,  KPTR,  KPUNCH 
GO  TO  (10,10,20,40)  10 

C*»»**»»INPUT 

10  READ  (KIN, 1000)  (X(I),  1=1, N) 

IF  (10  .EO.  1)  RETURN 

C»»«»«»»0UTPUT 

20  WRITE  (KOUT, 2000)  (X(I),  1=1, N) 

RETURN 

C«»«»»»»PUNCH 

40  WRITE  (KPUNCH, 3000)  (X(I),  1=1, N) 

RETURN 

1000  FORMAT  (8E10.0) 

2000  FORMAT  (1H  .1P10E13.3) 

3000  FORMAT  (1P8E10.3) 

END 


ooooo 
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SUBROUTINE  KVEClOdX, N,  10) 

BATCH  ORIENTED  INTEGER  VECTOR  I/O 
10=1  INPUT  ONLY 
10=2  INPUT  AND  OUTPUT 
10=3  OUTPUT  ONLY 
10=4  PUNCH 

DIMENSION  IX(1) 

COMMON  /INOU/  KIN,  KOUT,  KPTR,  KPUNCH 
GO  TO  (10,10,20,40)  10 

rt  tt  tt  ft  • ft  * iiip  or 

10  READ  (KIN.1000)  (IX(I) , 1=1 , N) 

IF  (10  .EQ.  1)  RETURN 

C*******OUTPUT 

20  WRITE  (KOUT, 2000)  (IX(I) , 1=1 , N) 

RETURN 

C^ftftftftft^pUliCH 

40  WRITE  (KPUNCH, 1000)  (IX(I),  1=1, N) 

RETURN 

1000  FORMAT  (8110) 

2000  FORMAT  (1H  ,10113) 

END 
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SUBROUTINE  LINEFD(KFIL,KOUNT)  ^ 

WRITES  KOUNT  LINEFEEDS  (SPACE  IN  COL  1)  ON  FILE  KFIL 

IF  ( KOUNT. LE.O)  RETURN 

DO  200  1=1, KOUNT 
WRITE  (KFIL, 1000) 

FORMAT  (1H  ) 

CONTINUE 

RETURN 

END 


- 110 
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C 

10 

1000 


SUBROUTINE  DAYTIM(KFIL) 

WRITES  THE  DATE  AND  THE  TIME  ON  FILE  KFIL 

CALL  TIME(LTIME) 

CALL  DATE(LDATE) 

WRITE  (KFIL, 1000)  LDATE,  LTIME 
FORMATUH  ,A10,2X,A10) 

RETURN 


WP 
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A1.  PROBLEM  DESCRIPTION 

GNPROP  is  an  interface  program  for  the  TIVAR  package  which  provides  a 
means  for  generating  time-varying  gains  by  backwards  integration  of  the 
Riccati  equation.  Although  GNPROP  runs  in  reverse-time,  those  gains  are 
punched  out  on  cards  by  GNPROP  in  forward  time,  and  may  be  inserted  directly 

into  the  TIVAR  input  deck.  In  this  appendix  the  procedure  for  running  the 
GNPROP  is  described. 
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A2.  TIME  VARIATIONS 

In  executing  GNPROP,  any  of  the  Input  parameters  can  be  changed  by  the 
user  at  any  time  t.  As  in  TIVAR,  there  are  two  methods  for  changing 
paraseters.  The  first  is  via  card  inputs;  the  second  is  via  a user  written 
subroutine  called  INTOLD.  (This  subroutine  corresponds  to  the  subroutine 
INTNEW  in  TIVAR,  but  because  GNPROP  operates  in  reverse-time,  it  is  named 
INTOLO.) 

A2. 1 System  Codes 

In  GNPP.CP,  a maximum  of  16  system  elements  may  be  changed  at  any  given 
time  step.  Each  element,  or  parameter,  is  identified  by  a unique 
alphanumeric  code  and/or  an  index  number.  Table  A1  defines  the  system  codes 
used.  As  in  TIVAR,  when  a given  parameter  I is  changed  at  time  t,  GNPROP 

sets  a flag  IFLAG(I)  equal  to  1 for  one  time  step.  The  implication  of  the 
parameter  change  is  then  addressed  via  the  internal  logic  in  GNPROP.  If 

there  have  been  no  parameter  changes,  the  internal  flags  remain  at  their 
nominal  zero  value . 

A2.2  Changing  System  Parameters 

Parameters  may  be  changed  at  time  t via  external  or  card  inputs.  The 
alphanumeric  code(s)  is  used  to  identify  the  parameter(s)  being  changed  at  a 
selected  time.  The  input  required  for  each  code  is  given  in  section  A3. 

Parameters  can  be  changed  periodically  via  ein  internal  subroutine 
INTOLD.  The  user  must  supply  his  own  code  — in  FORTRAN  IV  — to  use  this 
option.  For  any  specified  code  index  number,  I,  the  subroutine  INTOLD(KEY) 
is  called  once  every  NDT(I)  time  step  with  KEY=I.  The  16  values  of  NDT  are 
inputs  to  GNPROP.  When  INTOLD  is  called,  IFLAG(KEY)  is  set  to  1.  The  user 
must  supply  the  manner  in  which  the  parameters  are  to  be  changed.  If  no 
FORTRAN  code  is  supplied  to  update  the  parameters,  no  changes  are  made. 

Thus,  OIPROP  assumes  that  the  "new"  parameters  are  identical  to  the 
previously  existing  ones. 
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Table  A1:  PARAMETER  CODES  IN  GNPROP 


CODE 

INDEX 

DESCRIPTION 

A 

1 

System  A matrix 

B 

2 

System  B matrix 

C 

3 

Output  C matrix 

D 

4 

Output  D matrix 

E 

5 

Noise  matrix,  E 

QX 

6 

State  weightings  vector 

QY 

7 

Output  weightings  vector 

QU 

8 

Control  weightings  vector 

QR 

9 

Control  rate  weightings  vector 

RICI 

10 

First  guess  at  a riccati  solution 

QXI 

1 1 

State  weightings  vector  for  a discrete  time 

QYI 

12 

Output  wieghtings  vector  for  a discrete  time 

PRINT 

13 

Trigger  to  print  the  current  gains 

GNOUT 

14 

Trigger  to  punch  the  current  gains  ' 

DUMMY 

15-16 

Dummy  codes,  inactive  at  present 

- 115  - 
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A3.  INPUT  DECK  SETUP 

The  input  deck  structure  for  GNPROP  is  discussed,  along  with  the 
user-written  routine  INTOLD. 

A3. 1 Control  Cards 

There  are  five  major  control  cards  required  by  GNPROP. 

Card  1 - Title  Information 
Column  1:  blank 

Columns  2-80:  alphanumeric  title  information 

Card  2 - Gain  Output  Mode 

Card  3 - Dimension  Information,  515  Format 

Field  1 : NX  = number  of  system  states 

Field  2:  NX1=  number  of  noise  shaping  states 

Field  3:  NU  = number  of  control  Inputs  < 4 

Field  4:  NW  = number  of  random  noise  sources 

Field  5:  NY  = number  of  displayed  outputs 

The  size  restrictions  are; 


NX 


NU  < NDIM 
NY  < NDIM 


where  NDIM  is  the  array  size  in  the  DIMENSION  statemments  (15). 

Card  4 - Time  Information,  3E10.0  Format 

Field  1:  DEL  = discrete  time  step  interval  (sec) 

Field  2:  TO  = Initial  time  (sec) 

Field  3."  tend  = terminal  time  (sec) 

TO  and  TEND  are  specified  as  if  the  program  were  running  in 
forward  time.  They  correspond  to  the  same  quantites  in  TIVAR. 

Card  5 - NDT(I)  for  internal  time  breaks,  1615  Format 

The  16  fields  are  associated  with  the  16  parameter  codes  in  GNPROP  on  a 

one-to-one  basis.  The  I-th  field  is  associated  with  code  I.  NDT(I)  is  the 
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frequency  (number  of  time  steps)  at  which  subroutine  INTOLD  is  called 
internally  with  KEy=I,  starting  at  time  TEND.  Calling  INTOLD  with  KEY=I  sets 
flag  IFLAG(I)=1  for  one  time  step.  The  actual  parameter  value  must  be 
changed  by  user-written  code.  If  no  code  is  supplied,  the  associated 
parameters  retain  their  previous  values. 

A3. 2 External  Parameter  Cards 

The  remaining  input  cards  are  used  to  change  system  parameters  via 
external  read-in  at  specified  times.  The  deck  set-up  follows  the  same 
standard  form  as  for  TIVAR. 

Time  Card  - Cols.  1-4  Alphanumeric  TIME 

Cols.  11-20  Time  of  external  break,  E10.0  Format 
Code  Card  - Cols.  1-5  One  of  the  alphanumeric  codes  given  in 
Table  A1. 

Parameter  Cards  - The  new  parameter  values  required  by  the. 
associated  code. 

The  sequence  of  Code  Card  followed  by  new  parameter  values  is  repeated  for 
all  items  that  the  user  wishes  to  change  at  the  given  time.  To  change 
parameters  at  another  time,  input  a new  time  card  with  the  time  of  the 
desired  change,  followed  by  a code  card,  parameters  cards,  code  card,  etc. 
When  using  external  (card)  updates,  the  following  rules  must  be  observed: 

1.  Time  breaks  must  occur  in  decreasing  order,  since  the 
program  is  running  in  reverse-time . 

2.  The  code  cards  can  occur  in  any  order. 

3.  The  parameter  cards  must  immediately  follow  the  associated 
code  card. 

4.  Parameter  cards  must,  be  input,  as  the  program  expects 
to  read  in  new  values. 

5.  The  last  card  in  the  deck  is  an  end  card,  containing  the 
alphanumeric  END  in  cols.  1-3. 
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Table  A2  gives  the  required  input  for  each  of  the  active  codes,  as  well  as 
the  Initial  parameter  values  that  are  set  internally  by  the  program,  prior  to 
trTEND. 

Data  is  entered  on  the  cards  in  8E10.0  Format,  i.e.,  in  floating  point 
fields  of  10  columns  with  a maximum  of  8 fields  per  card.  The  numbers  may  be 

either  in  fixed-point  (decimal)  format  or  in  scientific  (exponential)  format  |' 

with  the  exponent  right-justified  in  the  field.  Matrices  are  entered  one  row  Lj 

r' 

at  a time . If  a row  contains  less  than  8 entries , the  remaining  fields  on  ! 

the  card  are  left  blank.  If  a row  contains  more  than  8 entries,  continue  on  i| 

I ' 

a second  card  for  that  row.  A new  row  always  begins  on  a new  card.  Vectors  i 

are  entered  in  similar  8E10.0  format:  the  first  entry  in  the  first  field,  the 
second  entry  in  the  second  field,  etc. 
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TABLE  A2:  CARD  DATA  INPUTS  AND  INITIALIZATION 


om. 

INDEX 

INPUT  DATA 

INITIAL  VALUUE 

A 

1 

*IJ  * 

1=1,..., NX, 

J=1,  .. 

. ,NX 

A=0 

B 

2 

®IJ  ’ 

1=1, ...,NX, 

j=1... 

. ,NU 

B=0 

C 

3 

"ij  i 

1=1,..., NY, 

j=1... 

. ,NX 

C=0 

D 

H 

“ij  ‘ 

1=1 NY, 

j=1... 

. ,NU 

D=0 

E 

5 

1=1, ...,NX, 

j=1... 

. ,NW 

E=0 

QX 

6 

Q 

xi  ’ 

1=1 NX 

Qx=0 

QY 

7 

^yi  • 

1=  1 , . . . , NY 

Qy  = 0 

QU 

8 

Q 

ui  ’ 

1=1, ...,NU 

Q =0 
u 

QR 

9 

Qri  ; 

1=1,..., NU 

Qp=0 

RICI 

10 

PINC 

* ^ s 1 

■*•>•••> 

NX+NU, 

j=l. 

. . . , NX+NU  PINC-0 

QXI 

11 

; 1=1, ...,NX 

o 

II 

X 

M 

cr 

QYI 

12 

; 1=1, . . . ,NY 

QI  =0 

PRINT 

13 

GNOUT 

14 

... 

T" 
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A4,  GNPROP  SUBROUTINES 

The  various  subroutines  that  constitute  the  GNPROP  package  are  discussed 
briefly,  along  with  the  named  common  blocks. 

A4. 1 Subroutine  Descriptions 

1.  Subroutine  GNPROP.  Provides  the  major  logic  control,  dimensioning 
declarations  and  Initialization  for  the  entire  package.  Virtually 
no  computation  is  done  in  GNPROP.  GNPROP  reads  the  5 control  cards. 

2.  Subroutine  DNDATE.  Reads  input  code  cards  and  the  subsequent 
data  input  cards.  Calls  INTOLD  periodically  as  per  NSTEP. 

Analagous  to  subroutine  UPDATE  in  the  TIVAR  program. 

3.  Subroutine  INTOLD.  User  written  routine  for.  internal  breaks. 

4.  Subroutine  GPFBGN.  Converts  Lc  and  outputs  equivalent 

T and  L . 

N opt 

5.  Subroutine  GNRITE.  Outputs  the  gains  as  specified  by  the  second 
control  card . 

A4.2  Common  Block  Usage 

Named  common  blocks  are  used  to  store  most  data  arrays  and  to  pass 
information  among  the  various  subroutines.  The  dimension  declarations  are 
given  in  the  primary  subprogram  GNPROP. 

1.  /OOMMUN/  NCODES,  IC0DES(16),  IFLAG(l6),  NSTEP(16) 

NCODES  = Number  of  possible  system  codes  in  Table  A1  (16) 

ICODES  = Alphanumeric  code  identifiers 

IFLAG  = 0 or  1 flags  to  indicate  pareuneter  changes 
NSTEP  = frequencies  for  internal  breaks,  control  card  5 

2.  /INOU/  KIN,  KOUT,  KPTR,  KPUNCH,  KDISK 

Logical  unit  numbers  for  input/output  devices 
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3.  /MAIN1/  NDIM,  ND1M1,  C0M1  /MAIN2/  COM2 

Common  blocks  required  for  library  subroutines.  NDIM  = dimension 
of  program  square  arrays. 

4.  /C0MP1/,  /C0MP3/ 

Common  blocks  used  for  internal  computations  and  internal  storage 
of  temporary  matrices. 

5.  /INPTS/  BD,  ADT 

Discretized  system  variables 

ADT  = discrete  system  augmented  A matrix  (transpose) 

BD  = discrete  system  B matrix 

6.  /INPTX/  N,  NX,  NX1,  NU,  AC 

Continuous  system  state  parameters 

AC  = augmented  system  A matrix  (codes  1 and  2) 

N = NX  + NU 

7.  /INPTX/  NY,  C /INPTW/  NW,  E 

Input  Information  for  display  outputs,  and  random  inputs 

8.  /TIMES/  TIME,  TGO,  TO,  TEND,  DEL 

TIME  = current  value  of  time,  t 

TGO  = time-to-go  = TEND  - TIME 

TO  = initial  time 

TEND  = final  time 

DEL  = discrete  time  step, 

9.  /WEIGHT/  QX,  QY,  QR,  P 

QX,  QY,  QR  = (augmented)  state,  output  and  control  rate  weightings 
respectively 

P = steady-state  Riccati  equation  solution. 

10.  /GAINBK/  CGD 

CGD  = state  feedback  gains 
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11.  /IHCRE/  QXI,  QYI,  PINC 

QXI,  QYI  = discrete  increments  to  the  cost  functional 

PINC  = increment  to  the  steady-state  Riccait  equation  solution 

12.  /STORE/  PTIME,  GSAVE 

Storage  for  the  time-varying  control  gains,  and  their  corresponding 
times. 
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A5.  SAMPLE  PROBLEM 


The  sample  problem  for  the  GNPROP  program  is  the  same  as  that  for  the 
TIVAR  program.  A description  of  this  problem  is  given  in  Section  6 of  this 
manual.  The  GNPROP  program  is  used  to  compute  the  (constant)  gains  for  the 
sample  problem.  This  section  contains  a listing  of  the  user  written 
subroutines,  the  input  data  deck,  and  a listing  of  the  output. 


A5. 1 User  Written  Subroutines  for  the  Sample  Problem 


C 

C 

C 


C 

C 


GNPP  - PROBLEM  DEPENDENT  SUBROUTINES  FOR  GPROPI 
INCLUDES: 

1 - SUBROUTINE  INTOLD 

SUBROUTINE  INTOLD(KEY) 

PERFORMS  INTERNAL  DOWNDATES 

PROBLEM  DEPENDENT  - TO  BE  SUPPLIED  BY  THE  USER 
COMMON 

/COMMUN/  NCODES,  IC0DES(l6),  IFLAG(l6),  NSTEP(16) 


IFLAG{.KBY)=1 

RETURN 

END 


! 
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A5-3  Output  Listing  for  the  Sample  Problem 

P-I-D  CONTROLLER.  GAIN  COMPUTATION 
5-Dec-76  18; 12 


DYNAMICS  READ  FROM  FILE:  PIDTI.DYN 

SQ-  TQtal  syst^  states  5 

NO.  Cff  NOISE  SHAPING  STATES  2 
NO.  OF  CONTROL  SYSTEM  INPUTS  1 
NO.  OF  RANDOM  NOISE  SOURCES  1 
NO.  OF  DISPLAYED  OUTPUTS  2 

INTEGRATION  TIME  STEP  = 0.050 

INITIAL  TIME  = -6.000 

TERMINAL  TIME  = -5.000 


INTERNAL 

DOWNDATES: 

INDEX 

CODE 

NDT 

13 

PRINT 

20 

14 

GNOUT 

20 

EXTERNAL 

DOWNDATE  AT 

TIME  = 

-5.000 

TGO  = 

0.000 

CODE:  QY 

QY  VICTOR: 

I.OOOE+06  0. 

EXTERNAL  DOWNDATE  AT  TIME  = 

OR  VECTOR; 

4. 22 IE-03 

DISCRETE  CONTROL  GAINS  AT  TIME 


-5 . 000  TGO  = 


0.000 


DGAIN  MATRIX: 
-1 .198E+01 
8.840E-»-00 


-2.904E+00 


-6.000 

4.102E+00  4.646E-01 


EQUIVALENT  CONTINUOUS  GAINS  LX,LU; 
CGAIN  MATRIX: 

-i:SSSi:81 


-3-335E-I-00  6.292E+00  6.467E-01 


CODE:  QR 


7.878E+00 


9.104E+00 


EQUIVALENT  CONTINUOUS  LOPT,  TN: 

L»,TN  MATRIX: 

-1.527E+00  -3.307E-01  6.239E-01  6.413E-02  9-027E-01 

9.916E-02 


i 
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The  following  is  a listing  of  the  gains  as  punched  by  the  GNPROP 
program.  These  cards  are  suitable  for  insertion  directly  in  the  TIVAR  input 
deck. 

time  -6.000 

1O8E-I-01  -2.904E+00  i*.102E+00  4.646E-01  7.878E+00 

8.840&I-00 


oooo oo  ooo  o o no oo 
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A6.  GNPROP  LISTING 
C NO  TABS,  COLONS,  OR  FORMFEEDS 

?0?£8Ses  •^^v^*^se-time  control  gain  propagation 

1 - block  data  gndat  - initializes  various  common  blocks 

2 - GNMAIN  - CALLS  SUBROUTINE  GNPROP 

3 - SUBROUTINE  GNPROP  - PRIMARY  SUBPROGRAM 


ALSO  REQUIRES  THE  FOLLOWING  SUBPROGRAM  FILES 

GNPP  - PROBLEM  DEPENDENT  SUBROUTINES  FOR  GNPROP 
INCLUDES 

1 - SUBROUTINE  INTOLD  - PERFORMS  INTERNAL  DOWNDATES 


1 


GNPIO  - I/O  FOR  GNPROP 
INCLUDES 

i : IBiSSHiai  dTO 

3 - SUBROUTINE  GNRITE 


SPECIFY  INTERNAL  DOWNDATES 
PERFORMS  EXTERNAL  DOWNDATES 
WRITES  OUT  THE  GAINS 


BLOCK  DATA  GNDAT 
INITIALIZES  VARIOUS  PARAMETERS 

COMMON 


DATA 

1 NCODES  /1 6/, 

2 ICODES 

2 /1HA,  1HB,  1HC,  1HD,  1HE,  2HQX,  2HQY,  2HQU,  2HQR,  4HRICI, 

2 3HQXI,  3HQYI,  5HPRINT,  5HGN0UT,  2»5HDUMMY/, 


^ ^I^,  KOU^,  KPTg,  KPUNC^,  KDIS|^ 


END 
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ssr: 


C 

c 


PROGRAM  GNMAIN 

1 (INPUT,  OUTPUT,  TAPE5=INPUT,  TAPE6=0UTPUT, 

2 PUNCH,  GNFILE,  TAPE7= PUNCH,  TAPE8=GNFILE) 

^LS^SU§I^TINE  GNPROP 

CALL  GNPROP 
END 
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SUBROUTINE  GNPROP 
PRIMARY  SUBPROGRAM 


1 

COMMON 

4 /MAIN1/  NDIM,  NDIM1 , COM1(15,15) 

5 /MAIN2/  C0M2(15,15) 

I 

8 /INPTS/  BD(15,4),  ADT(15,15) 

9 /INPTX/  N,  NX,  NX1,  NU,  AC(15,15) 

6 !• 

C /TIMES/  TIME.  TGO,  TO,  TEND,  DEL 
D /WEIGHT/  QX(30),  QY(30),  QR(30),  P(15,15) 

PINC(15,15) 

G /STORE/  PTIME(300),  GSAVE(15,300) 

DATA 

1 MNIL  /O/ 

SET  NDIM 

NDIMz15 

NDIM1=NDIM+1 

ZERO  THE  VECTORS  AND  MATRICES 

DO  10  1=1, NDIM 

QXlIUO.O 

QY(I)=0.0 

QR(I)=0.0 

8fiii);8:8 

DO  10  J=1,NDIM 
AC(I,J)=0.0 

?fy!;8:8 

PINC(I,J)=0.0 

CONTINUE 

NGAINrO 

SPECIFY  THE  TITLE 

READ  (KIN, 1120)  (TITLE(I) , 1=1,8) 

FORMAT  (8A1Q) 

IF  (EOF(KIN))  400,50 
CALL  PAGEFD(KPTR,1) 

WRITE  (KPTR,1125)  (TITLE(I),  1=1,8) 

FORMAT  (/.1H  .8A10) 

CALL  DAYTIM(KPTR) 

CALL  LINEFD(KPTR,2) 

SPECIFY  THE  GAINS  OUTPUT  MODE 
READ  (KIN .1160)  MODE 
WRITE  (KOOT,1135)  MODE 

FORMAT  (20H  GAINS  OUTPUT  MODE  =,‘^5,/,1H  ) 

GET  THE  PROBLEM  DIMENSIONS 

READ  (KIN, 1160)  NX,  NX1,  NU,  NW,  Ni 

FORMAT  (515) 

WRITE  (KPTR,1180)  NX.  NX1.  NU.  NW.  NY 
FORMAT  (28H  NO.  OF  TOTAL  SYSTEM  STATES  ,13,/, 
1 29H  no.  of  noise  SHAPING  STATES  ,12,/, 
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2 

i 

5 

C 

1200 


1225 

1 

2 

3 

C 


1310 


1320 


1330 

1 

65 

C 


C 

?0 

C 

80 


2090 


2092 


95 


2096 


29H  NO.  OF  CONTROL  SYSTEM  INPUTS, 12,/, 

29H  NO.  OF  RANDOM  NOISE  SOURCES  ,12,/, 

29H  NO.  OF  DISPLAYED  OUTPUTS  ,12,/, 

1H  ) 

SPECIFY  DEL,  TO,  AND  TEND 
READ  (KIN, 1200)  DEL,  TO,  TEND 
FORMAT  (3E10,0) 

TEND=T0+IFIX  ( ( TEND-T0-*-0 .001  )/DEL)  «DEL 
WRITE  (KPTR,1225)  DEL,  TO,  TEND 
FORMAT  (25H  INTEGRATION  TIME  STEP  = ,F10.3,/, 
17H  INITIAL  TIME  = ,F10.3,/, 

17H  TERMINAL  TIME  = !f10.3!/I 
1H  ) 


^^PFYtJN|ERNAL  downdates 


WRITE  (KPTR,1310) 

FORMAT  (22H  INTERNAL  DOWNDATES  ,19H  INDEX  CODE  NDT) 

DO  65  I=1,NCODES 

IF  (NSTEPU)  .EQ.  0)  GO  TO  65 

WRITE  (KPTR,1320)  I,  ICODES(I) , NSTEP(I) 

FORMAT  (26X,  12,  3X,  A5.  3X,  12) 

IF  (MDD(INT{(TEND-TO+0.6oi)/DEL),NSTEP(I))  .EQ.  0)  GO  TO  65 
WRITE  (KPTR,1330) 

FORMAT  (1H  ,16H  •••WARNING*** ,/ , 

1H  ,55H(TF-T0)  HAS  A NON-INTEGRAL  NUMBER  OF  INTERNAL  DOWNDATES) 
CONTINUE 

INITIALIZE  SOME  MORE  QUANTITIES 

NTERMS=5 

NXPUNX-i-1 

N=NX-*-NU 

TIMErTEND 

TGO=TEND-TIME 

TGNEXT=0.0 


MAIM  COMPUTATIONAL  LOOP  STARTS  HERE 

external  DOWNDATES 

SKIP  THIS  PART  IF  NO  PRINTING  WAS  CALLED  FOR,  OR  IF  TGO  IS 
IF  (IFLAG(13)  -EQ.  0)  GO  TO  100 
IF  (TGO  .LE.  0.0)  GO  TO  100 
CALL  TRANS2(N,NU,CGD,C0M2) 

WRITE  (KPTR,2090)  TIME 

FORMAT  (1H  .33HDISCRETE  CONTROL  GAINS  AT  TIME  = ,F10.3) 
CALL  MATlO(fc(5H2,NU,N,3) 

WRITE  (KPTR,2092) 

FORMAT  (IH  ,34HEQUIVALENT  CONTINUOUS  GAINS  LX,LU  ) 

CALE^EQUAti(QD,AC,N,N) 

CALL  MSCALE(QD(NXP1,1) ,C0M2,NU,N,-1.0) 

CALL  DSCRT(N,QD, DEL, COM2, COM 1,4) 

CALL  MS GALE (COM2, COM2, NU,N, DEL) 

CONTINUE 

CALL  MATIO(COM2.NU.N,3)  , , , , 

CALL  GM1NV(NU,n6,c6m^(1,NXP1) ,COM1(1 ,NXP1) ,MR,1) 

CALL  MMUL ( COM 1 ( 1 , NXP 1 ) , COM2 , NU , NU , NX , COM 1 ) 

WRITE  (KPTR,2096) 

CONTINUOUS  LOPT,  TN  ) 
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C SAVE  THE  GAINS 

100  IF  (IFLAG(14) .EQ.O  .AND.  TIME . GT .( TO+0 . 001 ) ) GO  TO  110 
IF  (TIME-0.001  .LE.  TO)  TIME=T0 
IF  (MODE  .EQ.  OhGO  TO  110 
IF  (TGO  .LE.  0.0)  GO  TO  110 
L=NGAIN«NU+1 
NGAIN=NGAIN+1 
PTIME(NGAIN)=TIME 
CALL  EQUATE(GSAVE(1 ,L) ,CGD,N,NU) 

IF  (L  .LT.  300)  GO  TO  110 
WRITE  (KPTR,2105) 

2105  FORMAT  MH  .37HINCREASE  STORAGE  ALLOCATION  FOR  GAINS) 
CALL  EXIT 


C 

110 


171 

172 


173 

C 

C 

174 


C 

175 


177 


178 


179 


C 


180 


PROPAGATE  THE  GAINS  FOR  ANOTHER  TIME  STEP 
IF  (TIME  ,EQ,  TO)  GO  TO  300 

IF  (IFLAGhO)  .EQ.  1)  CALL  MADD1  (N , N ,P, PINC,  P,  1 . 0) 

IF  (IFLAG(II)  .EQ.  0)  GO  TO  172 
DO  171  1=1, N 

IF  (IFLAG(12)  .EQ.  0)  GO  TO  174 
DO  173  1=1, NY 

COM1(J,I)=C1*C(I,J) 

CONTINUE 

CALL  MMULS(COM1,C,N,NY,N,P) 

DISCRETIZE  NEW  SYSTEM  DYNAMICS 

ADT,  BD  ARE  DISCRETIZED  AT  B OVER  THE  INTERVAL  (TIME-DEL, TIME) 
IF  (IFLAG(1)+IFLAG(2)  .EQ.  0)  GO  TO  175 
CALL  DSCRT(N,AC,DEL,ADT,C0M2,NTERMS) 

CALL  TRANS 1(N, ADT, ADT) 

CALL  EQUATE(BD,COM2(1 ,NXP1) ,N,NU) 

OBTAIN  DISCRETE  COST  WEIGHTS 
DO  177  1=1, NY 
C1=QY(I)»DEL/2.0 

5]*8(I,J) 

CALL  MAT2A(NY,N,C,QD,QD) 

DO  178  1=1. N 

QD(I,I)=QD(I,I)+QX(I)»DEL/2.0 

CONTINUE 

CALL  MMUL(ADT,QD,N,N,N,COM1) 

DO  179  1=1, NU 
C1=-1.0/(QR(T)*DEL) 

DO  179  J=1,N 

C0M2(J,I)=QD(J,I+NX)»DEL 

QXUD(J.I)=C1»C0M2(J,I) 

CONTINUE 

CALL  MAT6S(N,N,COM1,ADT,QD) 

CALL  MAT6S(N,NU,COM2,QXUD,QD) 


mihi 

CONTINUE 


UPDATE  CYCLE 

CALL  MMUL(P,BD,N,N,NU,COM1) 

CALL  MAT2A(N,NU,COM1,BD,CGD) 

DO  180  1=1. NU  , 

CGD(I.I)=CGD(I,I)+QR(I)*DEL 

CONTINUE 

CALL  GMINV(NU,NU,CGD,C0M2,MR,1) 
CALL  MMUL( COM 1. COM2, N.NU.NU.CGD) 
CALL  MAT5(CGD,6D,N,fcu|N,0oMl) 

DO  183  1=1, N 
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182 

183 


185 


C 

300 


C 

400 


DO  182  J=1,N 
COM1(I,J)=-C0M1(I,J) 

CONTINUE 

CALL  MMUL(C0M1,P,N,N,N,C0M2) 
DO  185  1=1, NU 


=C1»CGD(J,I) 


DO 
C^Ni 

PROPAGATE  CYaE 

TIME=TIME-DEL 

TGO=TGO^DEL 

CALL  EQUATE ( COM 1,ADT,N,N) 

CALL  MAT5S(QXUD,BD,N,NU,N,C0M1) 

CALL  EQUATE(P,QD,N,N) 

^L**^c2i&i?BxUD , QXUD , N , NU , - 1 . 0 ) 

CALL  MMULS(C0M1,CGD,N,N,NU,QXUD) 

CALL  EQUATE (CGD.QXUD.N.NU) 

GO  TO  70 

TIME  IS  EXPIRED  - WRITE  OUT  THE  GAINS  ON  A FILE 

NIrNU 

N2=N 

CALL  GNRITE(M0DE,NGAIN,N1,N2) 

DONE  - RETURN 

RETURN 

END 


onoo  ooooo 
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GNPIO  - I/O  FOR  GNPROP 
INCLUDES 

1 - SUBROUTINE  INTLET 


mmiM  giiii 


SUBROUTINE  INTLET 

SPECIFY  INTERNAL  DOWNDATES 

NSTEP(I)  REFERS  TO  THE  ITH  TYPE  OF  DOWNDATE 

NSTEP(I)=0  NO  INTERNAL  DOWNDATE 

NSTEP(I)=NDT  NDT  DELS  BETWEEN  SUCCESSIVE  INTERNAL  DOWNDATES 


1 

2 


COMMON 


/COMMUN/  NCODES.  ICODES(16).  IFLAG(16). 
/iNou/  KIN,  k6ut,  KPTR,  K^UNCH,  KDIS^ 


NSTEP(16) 


C READ  NSTEP(I)  FROM  ONE  CARD 

100  READ  (KIN, 1100)  (NSTEP(I),  1=1, NCODES) 

1100  FORMAT  (1615) 

RETURN 

END 
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C 


C 


100 

C 

110 

1010 

C 

120 

1030 

C 

130 

135 

C 

UO 

C 

150 

C 

160 

C 

1165 


SUBROUTINE  DNDATE(TGNEXT) 
PERFORMS  EXTERNAL  DOWNDATES 


1 5Mn/  NCODES.  IC0DES(16),  IFLAG(16),  NSTEP(16) 

2 /INOU/  KIN,  KOUT,  KPTR,  KPUNCH,  KDISK 
4 /MAIN1/  NDIM,  NDIM1 , C0M1(1) 

I ’• 

B /INPTtf/  NW 

C /TIMES/  TIME,  TOO,  TO,  TEND,  DEL 

DATA 

1 LEND,  LTIME,  LTGO  /3HEND,  4HTIME,  3HTGO/ 

NNX:NX*NDIM-«-1 

NXPUNX-fl 

TAKE  CARE  OF  INTERNAL  DOWNDATES 

DO  100  1=1, NCODES 

IFUG(I)=0 

IF  (NSTEpTiKEO.O)  go  TO  100  , ,, 

IF  (MOD(INT(TGO/DEL+0.01)  ,NSTEP(D)  .EQ.  0)  CALL  INTOLD(I) 
CONTINUE 

WRITE  (KPTR. 1010)  TGO 

FORMAT  (/,1H  ,23HEXT.  DOWNDATE  AT  TGO  = ,F10.3) 

SPECIFY  THE  TYPE  OF  EXTERNAL  DOWNDATE  (AND  TGNEXT  IF  LTGO) 
READ  (KIN, 1030)  IDEM,  BRKT 
FORMAT  (A5,5X,E10.0) 

IF  (EOF(KIN))  135,130 

CHECK  FOR  LEND 
IF  (IMIN.NE.LEND)  GO  TO  140 
TGNEXT=TEND-T0+0 . 1 
GO  TO  110 

CHECK  FOR  LTGO 
GO  TO  110 


i^^^folS^NkTJ^fME)  GO  TO  150 

TGNEXT=TEND-BRKT 

GO  TO  110 

DOME  WITH  SPECIAL  CHECKS 
CONTINUE 


SEARCH  THROUGH  THE  DOWNDATE  CODES,  ICODE(KEY) 
DO  160  KEY =1, NCODES 
IF  (IDEN.EQ.ICODES(KEY))  GO  TO  170 
CONTINUE 


CODE  WAS  ILLEGAL, 

WRITE  (KPTR, 1165)  IDEN 

FORMAT  (1H  ,5HCODE  ,A5,11H  IS  ILLEGAL) 

CALL  EXIT 
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SUBROUTINE  GNRITE(MODE,NGN,NU,N) 

WRITE  OUT  THE  (TRANSPOSED)  GAINS  IN  REVERSE  ORDER  FOR  TIVAR 

MODES  1 PUNCH 

M0DE=2  WRITE  ON  A FILE 


2 

4 

G 


1 


C 

10 


C 

100 


1020 

1040 

120 


COMMON 

/INOU/  KIN,  KOUT,  KPTR,  KPUNCH,  KDISK 
/MAIN1/  NDIM.  NDIM1,  COM1(1) 

/STORE/  PTliffi(300),  GSAVEU) 

DATA 

LTIME,  LDGAIN  /4HTIME,  5HDGAIN/ 


KsNGN-i-1 

L=NGN»NU«NDIM+1 

JEND=N«NDIM 

INITIALIZE  THE  FILE  UNIT  NUMBER  - KFIL 
IF  (MODE.LE.O)  RETURN 
IF  (M0DE.EQ.lj  KFILsKPUNCH 
IF  (M0DE.EQ.2)  KFIL=KDISK 
IF  (MDDE.GE.3)  RETURN 


OUTPUT.THE  GAINS 
L=L-NU«NDIM 


K=K-1 

IF  (K.EQ.O)  GO  TO  200 
TIME=PTIME(K) 

CALL  TRANSKn,NU,GSAVE(L)  ,COM1) 

WRITE  (KFIL,1020)  LTIME,  TIME,  LDGAIN 
FORMAT  (A4,6X,F10.3,/,A5) 

DO  120  1=1, NU 

WRITE  (KFiL104O)  (COMI(IJ),  U=I , JEND,NDIM) 
FORMAT  (1P8fel0.3) 

CONTINUE 
GO  TO  100 


C 

200 


DONE 

RETURN 
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LIBRARY  FDR  OPTIMAL  CONTROL  MODEL 


The  following  subroutines  from  the  Optimal  Control  Model  Library  are 
used  in  program  GNPROP,  as  well  as  in  TIVAR.  Listing  for  these  subprograms 
may  be  found  in  Section  7 of  this  manual. 


SUBROUTINE  DSCRT( N , A , DEL , EA , EAINT , NT) 
SUBROUTINE  GMINV(NR,NC,A,U,MR,MT) 

SUBROUTINE  MAT2(N1 ,N2,X. Y,Z) 
SUBROUTINE  MAT2A(N1 , N2, X, Y,Z) 
SUBROUTINE  MAT5(X,Y,N1 ,N2,N3,Z) 
SUBROUTINE  MAT5A(X, Y,N1 ,N2,N3,Z) 
SUBROUTINE  MAT6(N1,N2,X,Y,Z) 
SUBROUTINE  MMUL(X,Y,N1 ,N2,N3,Z) 

SUBROUTINE  MADD1(NR,NC,X,Y,Z,C1) 
SUBROUTINE  EQUATE(A,B,NR,NC) 
SUBROUTINE  MSCALE(A,B,NR,NC,C1 ) 
SUBROUTINE  TRANS1 (N ,A , AT) 

SUBROUTINE  TRANSP(N1 ,N2,X,XP0SE) 
FUNCTION  DOT (NR, A, B) 

FUNCTION  D0T2(NN,A,B) 

SUBROUTINE  VADD(N,C1,A,B) 

SUBROUTINE  VADD1(NN,C1 ,A,B) 
SUBROUTINE  VSCALE(X,Y,N,C1) 
SUBROUTINE  OUTPUT(M,N,X,U) 
SUBROUTINE  VECOUT(X,N,LET) 

SUBROUTINE  MATIO(X,NR,NC,IO) 
SUBROUTINE  VECTIO(X,N,IO) 

SUBROUTINE  KVECIO(IX,N,IO) 

SUBROUTINE  PAGEFD(KFIL.KOUNT) 
SUBROUTINE  LINEFD(KFIL,KOUNT) 
SUBROUTINE  DAYTIM(KFIL) 
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